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ABSTRACT 

Signals from radio pulsars show a wavelength-dependent delay due to dispersion in 
the interstellar plasma. At a typical observing wavelength, this delay can vary by tens 
of microseconds on five-year time scales, far in excess of signals of interest to pulsar 
timing arrays, such as that induced by a gravitational- wave background. Measurement 
of these delay variations is not only crucial for the detection of such signals, but 
also provides an unparallelled measurement of the turbulent interstellar plasma at au 
scales. 

In this paper we demonstrate that without consideration of wavelength- 
independent red-noise, 'simple' algorithms to correct for interstellar dispersion can 
attenuate signals of interest to pulsar timing arrays. We present a robust method for 
this correction, which we validate through simulations, and apply it to observations 
from the Parkes Pulsar Timing Array. Correction for dispersion variations comes at 
a cost of increased band-limited white noise. Wc discuss scheduling to minimise this 
additional noise, and factors, such as scintillation, that can exacerbate the problem. 

Comparison with scintillation measurements confirms previous results that the 
spectral exponent of electron density variations in the interstellar medium often ap- 
pears steeper than expected. We also find a discrete change in dispersion measure of 
PSR J1603-7202 of ~ 2 x lO'^ cm'^pc for about 250 days. We speculate that this has 
a similar origin to the 'extreme scattering events' seen in other sources. In addition, 
we find that four pulsars show a wavelength-dependent annual variation, indicating 
a persistent gradient of electron density on an au spatial scale, which has not been 
reported previously. 

Key vifords: pulsars: general — ISM: structure — methods: data analysis 



1 INTRODUCTION 

The fundamental datum of a pulsar timing experiment is 
the time of arrival (ToA) of a pulse at an observatory. In 
practise, the ToA is referred to the solar-system barycen- 
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tre in a standard time frame (e.g., barycentric coordinate 
time). This barycentric arrival time can be predicted us- 
ing a 'timing model' for the pulsar. The difference between 
the barycentric ToAs and the arrival times predicted by 
the timing model are termed residuals. The timing model 
can be refined using a least-squares fitting procedure to 
minimise the residuals, as performed by, e.g., the Tempo2 
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software (|Hobbs et al.|[200^ '). Since the timing model is al- 
ways incomplete at some level, we always see some level 
of post-fit residuals, which are typically a combination of 
'white' noise due to the uncertainty in the ToA measurement 
and 'red' (i.e., time-correlated) signal. For the majority of 
known pulsars the dominant red signal is caused by the in- 
trinsi c instability of th e pulsar, and termed 'timing noise' 
(e.g. , iHobbs et al.lboid ') . However, the subset of millisecond 
pulsars are stable enough that other r ed signals are poten- 
tially measurable jVerbiest et al.|[2009l ). Pulsar timing array 
projects, such as the P arkes Pulsar Timing Array (PPTA; 
[Manchester et al]l2012l ). aim to use millisecond pulsars to 
detect red signals su ch as: errors in the atomic time standard 



Hobbs et al.1l2012l): errors in the Solar System ephemeris 
201 0|); or the effect of Kravita.tional waves 



Champion et al.l— ..- ^ 

Yardlev et al.ll201(]| . l201ll : lvan Haasteren et al.ll201ll ). Each 



of these signals can be distinguished by the spatial corre- 
lation, i.e., how pulsars in different directions on the sky 
are affected. However, at typical observing wavelengths and 
time-spans, the variation of the dispersive delay due to tur- 
bulence in t he ionised interst ellar medium (ISM) dominates 
such signals (|You et al.ll2007l ). Fortunately for pulsar timing 
experiments, these delays can be measured and corrected 
using observations at multiple wavelengths. 
The dispersive group delay is given by 



teM = A 
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path 
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where A is the barycentric radio wavelengtlQ. The path in- 
tegral of electron density is the time-variable quantity. In 
pulsar experiments this is termed 'dispersion measure', DM, 
and given in units of cm~''pc. In principle, the instantaneous 
DM can be computed from the difference of two arrival times 
from simultaneous observations at different wavelengths, or 
more generally by fitting to any number of observations at 
more than one wavelength. 

The question of estimatio n and correction of DM(f) has 
previously been considered bv lYou et al] (|2007l ). They chose 
a 'best' pair of wavelengths from those available and esti- 
mated the DM at every group of observations. These ob- 
servation groups were selected by hand, as was the choice 
of wavelengths. Regardless of how the analysis is done, the 
estimated DM always contains white noise from differenc- 
ing two observations, and correcting the group delay al- 
ways adds that white noise to the arrival times. However 
the DM(t) variations are red, so they only need to be cor- 
rected at frequencies below the 'corner frequency' at which 
the power spectrum of the DM-caused fluctuations in group 
delay is equal to the power spectrum of the white noise in 
the DM(f) estimate. To minimise the additional white noise, 
they smoothed the DM(t) estimates over a time Ts to cre- 
ate a low-pass filter which cuts off the DM variations, and 
the associated white noise, at frequencies above the corner 
frequency. In this way, they avoided adding white noise at 
high frequencies where the DM-correction was unnecessary. 
Of course the added 'white' noise is no longer white; it is 
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white below the corner frequency, but zero above the corner 
frequency. 

Here we update this algorithm in two ways. We use 
all the observed wavelengths to estimate DM(f) and we in- 
tegrate the smoothing into the estimation algorithm auto- 
matically. Thus, the algorithm can easily be put in a data 
'pipeline'. We show the results of applying this new algo- 
rithm to the PPTA data set, whi ch is now about twice 
as long as when it was analysed bv lYou et al.l (|2007l V Ad- 
ditionally, we demonstrate that our algorithm is unbiased 
in the presence of wavelength-independent red signals, e.g., 
from timing noise, clock error, or gravitational waves; and 
we show that failure to include wavelength-independent red 
signals in the estimation algorithm will significantly reduce 
their estimated amplitude. 



2 THEORY OF DISPERSION REMOVAL 

We assume that an observed timing residual is given by 
toBS = icM + iDM(A/ARBF)^ where tcm is the common- 
mode, i.e., wavelength-independent delay and tDM is the dis- 
persive delay at some reference wavelength Aref. Then with 
observations at two wavelengths we can solve for both tcm 
and toM. 



turn = {toB. 

icM = (tOBS,2Ai 



toBS,2)AREF/(Al - A2), (2) 



toBS.iA2)/(A? -A^). (3) 



In a pulsar timing array, tcM would represent a signal of 
interest, such as a clock error, an ephemeris error, or the ef- 
fect of a gravitational wave. The dispersive component tuM 
would be of interest as a measure of the turbulence in the 
ISM, but is a noise component for other purposes. It is im- 
portant to note that toM is independent of tcM so one can 
estimate and correct for the effects of dispersion regardless 
of any common-mode signal present. In particular, common- 
mode red signals do not cause any error in torn- 

If more than two wavelengths are observed, solving for 
tcM and toM becomes a weighted least-squares problem, and 
the standard deviation of the independent white noise on 
each observation is needed to determine the weighting fac- 
tors. For wavelength i, we will denote the white noise by 
tw,i and its standard deviation by at so the observed timing 
residual is modelled as 



tOBS.i = tCM + tuui^i / ^REf)^ + tw,i- 



(4) 



The weighted least-squares solutions, which are minimum 
variance unbiased estimators, are 



toM = Aref ^ ^ l/crf ^ toBS.iAi /(Tj^ — 

i i 

^A?/a?^toBS,,/af)/A (5) 

i i 

= f a^/q-,^ y~^toBs.i/cr;^ — 

i i 

X^Ai/(7,^^toBS,iAi/cr,^^/A. (6) 

i i 

Here A is the determinant of the system of equations, 
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If one were to model only the dispersive term tuM, the 
weighted least-squares solution would become 



- _ -2 J2i toBS,iA?/o-,^ 

IDM — -^REF 



(7) 



However if a common- mode signal is present, this solution 
is biased. The expected value is 



2^1 \ I °"i 



(8) 



Some of the 'signal' torn is absorbed into tDM reducing the 
effective signal-to-noise ratio and degrading the estimate of 
DM. We will demonstrate this bias using simulations in Sec- 
tion g] 

It is important to note that the dispersion estimation 
and correction process is linear - the estimators t-ayi and 
fcM are linear combinations of the residuals. The corrected 
residuals toss, cor, i = ioss.i — (Ai/AREF)^teM, are also lin- 
ear combinations of the residuals. We can easily compute 
the white noise in any of these quantities from the white 
noise in the residuals. For example, we can collect terms 
in Equations ([5]) and ((G]) obtaining toM = X^i fiiioBS,i and 
icM = Yl,i bitoBS.i, where 

ai = ALF(A?/a?^l/a|-l/a?;^A,V'^|)/A (9) 

j j 

h = (l/a?^A>|-A?/a?^A,Va|)/A. (10) 

3 3 

Then, the white noise variances of the estimators can be 
written as CTtom = Y.i and a?tcM = Y.i blf^h 

The actual PPTA observations are not simultaneous at 
all frequencies, s o we cannot normally ap ply Equations ^ 
and ([6]) directly ^Manchester et al]|2012l ). We discuss how 
the least squares solutions for ioM and tcM can be obtained 
by including them in the timing model in the next section. 
However it is useful to have an analytical estimate of the 
power spectral density of the white noise that one can expect 
in these estimators and in the corrected residuals. At each 
wavelength Ai we have a series of Ni error estimates aij. The 
variance of the weighted mean is cr^i = 1/ X]j 1 / '^fj ■ This 
is the same as if we had a different number A'^ of observa- 
tions at this wavelength each of variance — o^iN . Thus, 
for planning purposes we can compute ami for each wave- 
length and conceptually resample each wavelength with an 
arbitrary number (A) of samples. Equations ([S]), ((6|, ((9]), 
and pO|) are invariant under scaling of all ai by the same 
factor so one can obtain the coefRcients ai and hi using ami 
in place of ai so the actual number (Ai) of samples need not 
enter the equations. 

If one had a series of A samples over a time span of 
ToBS each with variance a^ , the spectral density of the white 
noise would be = 2Tobs a'^ /N = 2roBS o"m- We can ex- 
tend this to a weighted white noise spectral density using 
the variance of the weighted mean. So the power spectral 
densities Pm,i play the same role as a'f in Equations ((5]), 
(O, (© and (jlOp . The coefficients {a^} and {hi} are func- 
tions of Ai and P^.i - Then we find Pu,,tdm ~ X^i o-iPw,i and 

Pw,TCM = X^i 6i Pw,i- 

Perhaps the most important property of these estima- 
tors is that Pw,TCM is less than or equal to the white noise 
spectrum of the corrected residuals Pw,cor.i in any band. 



Equality occurs when there are only two wavelengths. The 
values of P„,i, P„,cor,i, P™,tdm and P„,tcm are given for 
the PPTA pulsars in Table [T] Here Pu,,tdm is given at the 
reference wavelength of 20 cm. 

The situation is further complicated by red noise which 
depends on wavelength, but not as A^. For example, diffrac- 
tive angular scattering causes variations in the group de- 
lay, which scal e as the scat tered pulse width, i.e. approxi- 
mately as A* (|Rickettlll977l V Clearly such noise will enter 
the DM correction process. It can have the unfortunate ef- 
fect that scattering variations, which are stronger at long 
wavelengths, enter the short wavelength corrected residuals 
even though they are negligible in the original short wave- 
length data. This will be discussed in more detail in Section 

El 



3 DISPERSION CORRECTION TECHNIQUE 

Rather than solving for tcM and tuM for every group of 
observations, or re-sampling observations at each wave- 
length to a common rate, it is more practical to include 
parametrised functions for tcmit) and DM(t) in the tim- 
ing model used to obtain the timing residuals. To provide a 
simple and direct parametrisation we use piece-wise linear 
models defined by fixed samples tcM{tj) and DM(tj) for j = 
1,...,A,. 

It is also required to introduce some constraints into the 
least-squares fitting to prevent covariance with other model 
parameters. For example, the values of DM{tj) are natu- 
rally covariant with the mean dispersion measure parame- 
ter, DMo, which is central to the timing model. To eliminate 
this covariance, we implement the linear equality constraint 
that J2i^i'DM{tj) = 0. Additionally, the series tcM(^j) is 
covariant with the entire timing model, however in practise 
the sampling interval is such that it responds very little to 
any orbital parameters (in the case of binary systems). We 
constrain tcM(ij) to have no response to a quadratic poly- 
nomial, or to position, proper motion, and parallax. These 
constraints are implemented as part of the least-squares fit 
in Tempo2, as described in Appendix 1X1 

The ch oice of sampling interval, Tg is essentially the 
same as in IYou et al] (120071 ). The process of fitting to a 
piece-wise linear function is equivalent to smoothing the 
T)M{t) time series with a triangle function of base 2Ts. 
This is a low pass filter with transfer function -fftri(/) = 
(sin(7r/rs)/7r/rs)^. We adjust Ts such that the pass band 
approximately corresponds to the corner frequency fc at 
which the power spectrum of the DM delays, Ptdm, ex- 
ceeds that of the white noise, Piu.tdm- Note that this corner 
frequency is independent of reference wavelength at which 
tDM is defined. 

To determine this corner frequency we need an estimate 
of the power spectrum of toM, so the process is inherently 
iterative. We can obtain a first estimate of Ptdm(/) from 
the diffractive time scale, rdift , at the reference wavelength. 
For signals in the regime of strong scattering, which includes 
all PPTA observations, Tdiff is the time scale of the diffrac- 
tive intensity scintillations. For the PPTA pulsars, Tdiff is 
usually of the order of minutes and can be estimated from 
a d ynamic spectrum ta ken during normal observations (see 
e.g. ICordes et al11l990l ). 
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Table 1. The estimated power spectral density before {Pw) and after (P^j.cor) correction of the white noise for each PPTA pulsar at 
each of the three wavelengths, and the expected white noise power spectral density in the 'common mode' signal (P^^tcm) ^'^d in (dm 
at 20 cm (Pii,,tdm)i all expressed relative to the power spectral density of the uncorrected 20-cm residuals. Also shown is the effect of 
optimising the observing time, expressed as the ratio of Pu.'.TCM estimated for optimal observing and P^j^tcm with the cmrent observing 
strategy (q = 0.5), and Qopt the optimal fraction of time spent using the dual 10- and 50cm-cm observing system. 
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Rather than directly compute Ptdm, it is attractive to 
begin with the structure function, which is a more conve- 
nient statistic for turbulent scattering processes and is more 
stable when only a short duration is available. The structure 
function of tuM is given by 

DtDmM = {{tDM{t) ~ toM{t + T)f) = {\/2ncfD^{T), 

(11) 

where D^(t) is the phase structure function. If we assume 
that the electron density power spectrum has an expo- 
nent of -11/ 3, i.e., Kolmogorov tu rbulence, then D^{t) = 
(r/rdiff)^''3 (|Foster fc Cordeslll990l ). The structure function 
Dt-dm{t) can therefore be estimated from Tdiff, or directly 
from the t-Dviit) once known. 

As described in Appendix |B] we can use the structure 
function at any time lag r to obtain a model power spectrum 
using 

PtDm(/) ^ 0.0112 I>TDM(r)r-''^''3(gpy-,-l/3^-8/3 (^2) 

The term (spy) is the number of seconds per year. Here 
Dtum is in s2, r in s, / is in yr^^ and Ptdm is in yr^. 

The spectrum of the white noise can be estimated from 
the ToA measurement uncertainties as discussed in section 
2. However, often there are contributions to the white noise 
that are not reflected in the measurement uncertainties and 
so we prefer to estimate directly from the power spec- 
trum of the residuals. 



4 TEST ON SIMULATED OBSERVATIONS 

When dealing with real data sets it is not trivial to show 
that the DM-corrected residuals are 'improved' over sim- 
ply taking residuals from the best wavelength (|You et al.l 



[20o3). This is because much of the variations in DM are ab- 
sorbed into the fit for the pulsar period and period deriva- 
tive. Therefore the root-mean-square deviation (RMS) of the 
residuals from a single wavelength may not decrease signif- 
icantly even though the RMS of the DM(f) variations that 
were removed is large. To demonstrate that the proposed 
procedure can estimate and remove the dispersion, and that 
it is necessary to include the common-mode in the process, 
we perform two sets of simulations. 

The observing parameters, i.e., Tobs, Ni, Oij, -Ddm(t), 
of both simulations are based on the ob servations of PSR 
J1909 -3744 in the PPTA 'DRl' data set (|Manchester et al.1 
l2012h . We find it useful to demonstrate the performance of 
the DM correction process in the frequency domain, but it 
is difficult to estimate power spectra of red processes if they 
are irregularly sampled. Therefore we first use simulations 
of regularly sampled observations with observing parame- 
ters similar to those of PSR J1909— 3744 to demonstrate 
the performance of the DM correction algorithm. Then we 
will simulate the actual irregularly sampled observations of 
PSR J1909— 3744 to show that the ultimate performance of 
the algorithm is the same as in the regularly sampled case. 

4.1 Regular sampling, equal errors 

We will compare the power spectra produced after fitting 
for DM(t) with and without simultaneously fitting for a 
common-mode signal. To generate the simulated data sets, 
we first generate idealised ToAs that have zero residual from 
the given timing model. Then we add zero-mean stochastic 
perturbations to the ideal ToAs to simulate the three compo- 
nents of the model: (1) independent white noise, correspond- 
ing to measurement error; (2) wavelength independent red 
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noise, corresponding to the common-mode; (3) wavelength 
dependent red noise representing DM(t). 

We simulate the measurement uncertainty with a white 
Gaussian process, chosen to match the high frequency power 
spectral density of the observed residuals. The simulated Pw 
is 2.2 X 10"^^ 4.3 X 10~^° and 2.6 x 10"^'' yr^ at 10, 20 and 
50cm cm respectively. For the common mode we choose a 
Gaussian process with a spectrum chosen to match a com- 
mon model of the incoherent gravitational wave background 
(GWB), i.e. Pawsd) = (^^wB/lSTr')/-"''^' (Ijenet et all 
l2006l : lHobbs et al.ll2009l '). For the DM we use a Gaussian pro- 
cess with a power spectrum Pdm(/) = ^dm/~*^^, where 
AuM is chosen to match the observed DM fluctuations in 
PSR J1909-3744 shown in Figures [S] and [51 and the spectral 
exponent is chosen to match that expected for Kolmogorov 
turbulence (|Foster fc Cordes|[l990l ). The levels of Ptdm and 
PawB are similar so that the same sample intervals can be 
used for both DM(ti) and tcM{ti), but this is not necessary 
in general and will not always be desirable. 

For both algorithms, we estimate the pre- and post- 
correction power spectra of the 20-cm residuals in four noise 
regimes: P„; P^u+Pdm; Pw + Pgwb; and Pu, -|-Pdm + Pgwb. 
In order to minimise the statistical estimation error, we aver- 
age together 1000 independent realisations of the spectra for 
each algorithm. We note that although the averaged power 
spectra suggest that the input red noise signals are large, 
the noise on a single realisation is such that the red signals 
are at the limit of detection. To illustrate this, the 90% con- 
fidence limits for both the 1000 spectrum average and for a 
single realisation, are shown on the power spectra in Figures 
[T]and[21 

We show the effect of using the interpolated model for 
DM(t), but not fitting for the common-mode signal tcM{t), 
in Figure [1] This algorithm is well behaved when the GWB 
is not present, as shown in the two lower panels. In this 
case the DM correction algorithm removes the effect of the 
DM variations if they are present and increases the white 
noise below the corner frequency by the expected amount. 
Importantly, when the model GWB is included, i.e., in the 
two top panels, a significant amount of the low-frequency 
GWB spectrum is absorbed into the DM correction. This is 
independent of whether or not DM variations are actually 
present because the DM correction process is linear. 

We show the full algorithm developed for this paper, 
using interpolated models for both DM(t) and the common- 
mode signal tcM{t), in Figure O One can see that the algo- 
rithm removes the DM if it is present, regardless of whether 
the GWB is present. It does not remove any part of the 
GWB spectrum. When the GWB is not present, as shown 
in the two lower panels, the algorithm remains well behaved. 
As expected, it increases the white noise below the corner 
frequency by a larger factor than in the previous case. This 
is the 'cost' of not absorbing some of the GWB signal into 
the DM correction. Although it has a higher variance than 
for the previous case, our DM{t) is the lowest variance un- 
biased estimator of the DM variations in the presence of 
wavelength-independent red noise. This increase in white 
noise is unavoidable if we are to retain the signal from a 
GWB, or indeed any of the signals of interest in PTAs. 

The power spectra presented in Figure [5] demonstrate 
that the algorithm is working as expected, in particular that 
it does not remove power from any wavelength-independent 



signals present in the data. We note, however, two limita- 
tions in these simulations: the regular sampling and equal 
errors are not typical of observations, nor have we shown 
that the wavelength-independent signal in the post-fit data 
is correlated with the input signal (since our power spec- 
trum technique discards any phase information) . These lim- 
itations will be addressed in the next section. 



4.2 Irregular sampling, Variable error bars 

In order to test the algorithm in the case of realistic sampling 
and error bars, we repeated the simulations using the actual 
sampling and error bars for pulsar J1909— 3744 from the 
PPTA. We use the same simulated spectral levels for the 
GW and DM as in the previous section. The results are also 
an average of 1000 realisations. 

As a direct measure of performance in the estimating 
DM(t), we compute the difference between the DM esti- 
mated from the fit to the residuals, DMcst(i), and the DM 
input in the simulation, DMin(t). To better compare with 
the timing residuals, we convert this error in the DM into 
the error in toM{t) at 20cm using Equation ([T|). Note that, 
although the residuals were sampled irregularly, the original 
DMin(t) was sampled uniformly on a much finer grid. Fur- 
thermore, the estimated DMcst (t) is a well defined function 
that can also be sampled uniformly. Thus it is easy to com- 
pute the average power spectrum of this error in tDM(f) as is 
shown in Figure [3] We also plot the spectrum of the initial 
white noise, and the spectrum of the white noise after cor- 
rection. If the algorithm is working correctly the white noise 
after correction should exactly equal the error in toM{t) plus 
the white noise before correction, so we have over plotted the 
sum of these spectra and find that they are identical. 

The spectrum of the error in toMit) shows the expected 
behaviour below the corner frequency. Above the corner fre- 
quency (where the correction is zero), it falls exactly like the 
spectrum of tDM(^) itself, i.e., as f~^^^. By comparing the 
right and left panels one can see that the DM correction is 
independent of the GWB. 

We can also demonstrate that the model GWB sig- 
nal is preserved after DM correction, by cross-correlation 
of the input model GWB with the post-correction residu- 
als. If the GWB signal is preserved this cross-correlation 
should equal the auto-correlation of the input GWB signal. 
We show the auto-correlation of the input and four different 
cases of the cross-correlation of the output in Figure 2] The 
cross-correlations are for two bands (20 and 50cm cm shown 
solid and dashed respectively), and for two different fitting 
algorithms (with and without tcM (t) shown heavy and light 
respectively). Again it can be seen that, without fitting for 
the common-mode tcM{t), a significant portion of the GWB 
is lost. In fact, it is apparent from the large negative correla- 
tion at 50cm cm that the 'lost' power is actually transferred 
from the 20-cm residuals to those at 50cm cm. Although it 
may be possible to recover this power post-fit, it is not clear 
how to do this when the GWB and DM signals are unknown. 
Finally, we note that when the common mode is used, the 
50cm-cm residuals preserve the GWB just as well as the 20- 
cm residuals, even though they carry the majority of the 
DM(t) variation. 



6 M. J. Keith et al. 




Frequency (yr ) 



1 10 
Frequency (yr ^) 



Figure 1. Average power spectra of pre- and post-correction timing residuals, in ttie 20-cm band, with four combinations of signals. 
The solid line shows the pre-correction spectrum and dashed line shows the post-correction spectrum. For the cases where variations 
in DM are included in the simulation, the pre-correction spectrum without DM variations is shown with a dotted line. Here the fitting 
routine uses the DM(t) interpolated fitting routine, without fitting a common-mode signal. The vertical bars on the left of each panel 
show the 90% spectral estimation uncertainty for a single realisation (left- most bar) and the average of fOOO realisations (right bar). 




Lag (days) 

Figure 4. Cross-correlation of post-correction residuals with 
input model GWB, for the simulations representing PSR 
J1909— 3744. Solid lines show data from the 20-cm wavelength 
and dashed lines show data from the 50cm-cm band. The cor- 
rection was computed with and without fitting for the common 
mode, indicated by thick and thin lines respectively. The auto- 
correlation of the input GWB is plotted as a dotted line, but it 
is completely obscured by the heavy solid line for the cross cor- 
relation in the 20-cm band. 



4.3 The Robustness of the Estimator 

The proposed DM correction process is only optimal if the 
assumptions made in the analysis are satisfied. The pri- 
mary assumptions are: (1) that there is an unmodelled 
common-mode signal in the data; (2) that the residuals 
can be modelled as a set of samples toi(ij) = tcmitj) + 
tDM(ti)(Ai/AREF)^ (3) the variances of the sam- 

ples twiitj) are known. 

If Assumption 1 does not hold and we fit for tcyi{tj), 
then our method would be sub-optimal. However in any pul- 
sar timing experiment, we must first assume that there is 
a common-mode signal present. If tcM(ij) is weak or non- 
existent then we will have a very low corner frequency and 
effectively we will not fit for tcmitj)- So this assumption is 
tested in every case. 

Assumption 2 will fail if there are wavelength depen- 
dent terms which do not behave like A^, for example the 
scattering effect which behaves more like A^. If these terms 
are present they will corrupt the DM estimate and some 
scattering effects from longer wavelengths may leak into the 
shorter wavelengths due to the correction process. However 
the correction process will not remove any common-mode 
signal, so signals of interest to PTAs will survive the DM 
correction unchanged. 

Assumption 3 will not always be true a priori. Recent 
analysis of single pulses from bright MSPs have shown that 
pulse-to-pulse variations contribute significant white noise in 
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Figure 2. As for Figure[T] except the fitting routine uses tlie DM(t) interpolated fitting routine in addition to the wavelength independent 
signal, C(t). 
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Figure 3. Average power spectra of the error in DM(t) after fitting to simulations with realistic sampling and uncertainties. The 
simulations contained white noise, DM variations and, in the left panel, a model GWB. The solid black line shows the power spectrum 
of DMeat(i) — DMin(i). The dotted line is the power spectrum of the white noise only. The dashed line is the post-correction power 
spectrum of the residuals, after subtracting the model GWB signal if present. The crosses mark the sum of the black line and the dotted 
line. 
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Table 2. Scintillation and dispersion properties for the 20 PPTA 
pulsars, at a reference wavelength of 20 cm. The scintillation 
bandwidth (I'o) and time scale {r^iff) are averaged over a large 
number of PPTA observ ations except for values in parenthesis 
which are are taken from I You et al.l ll2007t l. Diooo is the value of 
the structure function at 1000 days and Tg is the optimal sampling 
interval for tDM(i). 



Source 


1^0 


Tdifl 


£'1000 






(MHz) 


(s) 







J0437-4715 


1000 


2486 


1.6 


0.2 


J0613-0200 


1.64 


4500 


0.3 


1 


J0711-6830 


36 


1962 


1.9 


2 


J1022+1001 


65 


2334 


0.14 


2 


J1024-0719 


268 


4180 


6.2 


1 


J1045-4509 


(0.094) 


(119) 


690 


0.25 


J1600-3053 


0.09 


271 


24 


0.5 


J1603-7202 


5 


582 


5.5 


1 


J1643-1224 


0.022 


582 


65 


0.5 


J1713+0747 


24 


2855 


0.31 


1 


J1730-2304 


12.4 


1615 


20 


1 


J1732-5049 


5.4 


1200 


10.0 


1 


J1744-1134 


60 


2070 


1.3 


1 


J1824-2452A 


(0.025) 


(75) 


250 


0.33 


J1857+0943 


5.5 


1464 


0.9 


2 


J1909-3744 


37 


2258 


3.5 


0.33 


J1939+2134 


1.2 


327 


8.9 


0.33 


J2124-3358 


(1170) 


(10705) 


0.4 


2 


J2129-5721 


17.1 


3060 


0.49 


2 


J2145-0750 


195 


3397 


0.15 


2 



excess of th at expected from the formal ToA measuremen t 
uncertainty jOslowski et al.ll201ll : rShannon fc Cordes|[20l3 ) . 
Indeed, many pulsars appear to show some form of addi- 
tional white noise which is currently unexplained but could 
be caused by the pulsar , the interstellar medium, or the ob- 
servi ng system (see e.g. ICordes fc Down3ll985l : iHotan et al.l 
I2OO5I ). In any case, we cannot safely assume that the un- 
certainties of the timing residuals aij accurately reflect the 
white noise level. If the aij are incorrect, our fit parame- 
ters, tc{t) and td{t), will no longer be minimum variance 
estimators; however, they will remain unbiased. This means 
that DM estimation will be unbiased and the DM correction 
will not remove any GWB (or other common-mode signal) 
although the correction process may add more white noise 
than optimal. It should be noted that if all the aij were 
changed by the same factor our DM correction would be 
unchanged. Fortunately the actual white noise is relatively 
easy to estimate from the observations because there are 
more degrees of freedom in the white noise than in the red 
noise, so in practise we use the estimated white noise rather 
than the formal measurement uncertainties an. 



5 APPLICATION TO PPTA OBSERVATIONS 

We have applie d the new DM correct ion technique to the 
PPTA data set (jManchester et al.l20l3 ). Observations of the 
PPTA are made in three wavelength bands: '10 cm' (~ 3100 
MHz); '20cm' (~ 1400 MHz); and '50cm cm' (~ 700 MHz). 
The 10-cm and 20-cm bands have been constant over the 
entire time span, however the long wavelength receiver was 
switched from a centre frequency of 685 MHz to 732 MHz 



around MJD 55030 to avoid RFI associated with digital TV 
transmissions. To allow for changes in the intrinsic pulse 
profile between these different wavelength bands, we fit for 
two arbitrary delays between one wavelength band and each 
of the other bands. However we did not allow an arbitrary 
delay between 685 and 732 MHz because the pulse shape 
does not change significantly in that range. 

We began our analysis by using the procedure described 
in Section [3] to compute pilot estimations of DM(t) and 
tcm{t) for each of the 20 pulsars, using a sampling inter- 
val Ts = 0.25 yr. Figure [5] shows the DM(i) derived from the 
above. Our results are c onsistent with the measurements 
made bv lYou et af] (|2007l ) for the ~ 500 days of overlapping 
data, which is expected since they are derived from the same 
observations. 



5.1 Determining the sampling interval 

As discussed in Section [3l we can use the diffractive time 
scale Tdiff to predict the magnitude of the DM variations in 
a given pulsar. This value can be computed directly from 
observations, however it is always quite variable on a day to 
day time scale (see Section [6] for discussion), and for a few 
pulsars Tdiff approaches the duration of the observations, so 
it can be hard to measure. Nevertheless we have obtained 
an estimate of the average Tdiff from the dynamic spectra 
for each pulsar, and this is given in Table [2l We have not 
provided an error estimate on the average Tdiff because the 
variation usually exceeds a factor of two so the values tabu- 
lated are very rough estimates. 

We also computed the structure function Dtdm directly 
from the tDM values. These structure functions, scaled to 
delay in /is^ at 20 cm, and those estimated from Tdiff, are 
shown in Figure [S] The value Dtdm (1000 days) is given in 
Table [1 

For each pulsar we also make an estimate of the white 
noise power directly from the power spectrum of the resid- 
uals. The estimates of at each wavelength are given in 
Tabled 

We then use the Dtdm estimates and Equation (|12|l to 
generate a model power spectrum Ptdm(/) at the reference 
wavelength (20 cm) for each pulsar. These assume a Kol- 
mogorov spectral exponent. From these model spectra and 
the corresponding Pw,tdm, tabulated in Table [1] we deter- 
mine the corner frequency and the corresponding sample 
interval Ts for DM for each pulsar. As we do not have any 
a priori knowledge of tcm for the PPTA pulsars we choose 
the same sample interval for tcM as for toM- 

5.2 Results 

The measured DM(t) sampled at the optimal interval Ts 
is overlaid on the plot of the pilot analysis with Ts — 
0.25 yr on Figure [S] It is not clear that there are measur- 
able variations in DM in PSRs J1022-I-1001, J2124-3358 or 
J2145— 0750, but one can see that there are statistically sig- 
nificant changes with time for the other pulsars. In general, 
the 'optimally sampled' time series (dashed line) follows the 
DM trend with less scatter. However, there are some sig- 
nificant DM fluctuations that are not well modelled by the 
smoother time-series. In particular we do not model the sig- 
nificant annual variations observed in PSR J0613— 0200, and 
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Figure 5. DM as a function of time for 20 PPTA pulsars. Solid lines show values measured with intervals of 0.25 yr. In the cases where 
the optimal Ts at a wavelength of 20cm is longer than 0.25 yr, a dashed line is added showing DM(ti) measured with this time step. 
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Figure 6. Structure functions of dispersive delay at 20 cm. The square markers indicate the structure function as measured directly 
from the DM time-series in Figure |5] error bars are derived by simulation of white noise. The solid lines show the extrapolation from 
the scintillation time scale Tdiff assuming a Kolmogorov spectrum, dashed lines mark the region occupied by 68% of simulated data sets 
having Kolmogorov noise with the same amplitude. These lines indicate the uncertanty in the measured structure functions resulting 
from the finite length of the data sets. The dotted lines show a Kolmogorov spectrum with the amplitude set to match the real data at 
a lag of 1000 days. 
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Table 3. Impact of the DM corrections on the timing parame- 
ters, as determined in the 20-cm band. For each pulsar we present 
the change in u and u due to the DM correction, relative to the 
measurement uncertainty, and the ratio of the RMS of the residu- 
als before (Eprc) and after (Spost) DM correction. Also included 
is the ratio of the power spectral density before (Ppre) and after 
(^post) correction, averaged below fc- The final column indicates 
if we believe that the DM corrections have 'improved' the data 
set for the purpose of detecting common- mode red signals. 



PSR 




|Ai> 


Spost 


-Ppost 


Imp. 








Ppro 


J0437-4715 


92 


48 


0.6 


0.15-0.25 


Y 


J0613-0200 


0.16 


2.9 


1.1 


0.3-1.2 


y 


J0711-6830 


3.9 


5.5 


1.0 


0.4-1.6 


y 


J1022-I-1001 


1.4 


0.3 


1.0 


0.6-2.6 


n 


J1024-0719 


1 


0.91 


1.0 


0.2-0.7 


Y 


J1045-4509 


28 


11 


0.7 


0.22-0.39 


Y 


J1600-3053 


35 


0.51 


1.0 


0.4-0.8 


Y 


J1603-7202 


2.4 


2.5 


1.0 


0.2-0.9 


Y 


J1643-1224 


11 


0.73 


1.7 


1.3-3.1 


N 


J1713-I-0747 


3.2 


6.2 


1.0 


0.2-0.7 


Y 


J1730-2304 


6.5 


1.8 


1.1 


0.9-3.2 


n 


J1732-5049 


2.6 


2.8 


1.0 


0.4-1.4 


y 


J1744-1134 


5.4 


0.48 


1.0 


0.5-2.0 


n 


J1824-2452A 


24 


31 


0.7 


0.29-0.56 


Y 


J1857-I-0943 


4.3 


1 


1.0 


0.2-1.0 


y 


J1909-3744 


28 


5 


1.0 


0.44-0.79 


Y 


J1939-I-2134 


13 


1.7 


0.7 


0.34-0.67 


Y 


J2124-3358 


0.25 


0.056 


1.0 


0.5-1.9 


y 


J2129-5721 


3 


2.1 


1.1 


0.7-2.8 


n 


J2145-0750 


0.22 


0.18 


1.0 


0.2-1.0 


y 



we must add a step change to account for the 250 day in- 
creeise in DM observed in PSR J1603— 7202 (these features 
are discussed more fully in Section [8]). These variations do 
not follow the Kolmogorov model that was used to derive the 
optimal sampling rate, and therefore we must use a shorter 
Ts so we can track these rapid variations. These results il- 
lustrate the importance of making a pilot analysis before 
deciding on the sample interval. The ISM is an inhomoge- 
neous turbulent process and an individual realisation may 
not behave much like the statistical average. The DM(t) for 
PSR J1909— 3744 is also instructive. It is remarkably linear 
over the entire observation interval. This linearity would not 
be reflected in the timing residuals at a single wavelength 
because a quadratic polynomial is removed in fitting the 
timing model. It can only be seen by comparing the residu- 
als at different wavelengths. Such linear behaviour implies a 
quadratic structure function and a power spectrum steeper 
than Kolmogorov. 

5.3 Performance of DM Correction 

The simplest and most widely used metric for the quality of 
timing residuals is the RMS of the residuals. Thus a natural 
measure of the performance of DM correction would be the 
ratio of the RMS of the 20-cm residuals before and after 
DM correction. This ratio is provided in Table [S] However, 
for most of these pulsars, the RMS is dominated by the 
white noise and so does not change appreciably after DM 
correction. Furthermore much of the effect of DM(t) varia- 



tions is absorbed by fitting for the pulse frequency and its 
derivative. Thus the ratio of the RMS before and after DM 
correctio n is not a very se nsitive performance measure. As 
noted bv lYou et al.] (|2007h . the DM correction has a signifi- 
cant effect on the pulsar spin parameters, which can give an 
indication of the magnitude of the DM correction. Table [3] 
lists the change in v and as a factor of the measurement 
uncertainty, caused by applying the DM correction. How- 
ever, there are systematic uncertainties in the estimation of 
the intrinsic values of v and v that may be greater than the 
error induced by DM variations. 

Judging the significance of the DM corrections depends 
on the intended use of the data set. Since a major goal of the 
PPTA is to search for common-mode red signals, we choose 
to consider the impact of the DM corrections on the low fre- 
quency noise. In principal, the DM correction should reduce 
the noise at frequencies below /c, and therefore we have 
estimated the ratio of the pre- and post-correction power 
spectrum of the 20-cm residuals, averaged over all frequen- 
cies below /c. We caution that the spectral estimates are 
highly uncertain, and for many pulsars we average very few 
spectral channels so the error is non-Gaussian. Therefore, 
we present these ratios in Table [3] as an estimated 68% un- 
certainty range, determined assuming the spectral estimates 
are x^-distributed with mean and variance equal to the mea- 
sured mean power spectral density. 

There are 9 pulsars for which the DM correction ap- 
pears to significantly reduce the low frequency noise, and 
therefore increases the signal-to-noise ratio for any common- 
mode signal in the data. These pulsars are listed with a 'Y' 
in Table [31 There are 10 pulsars for which the change in low 
frequency power is smaller than the uncertainty in the spec- 
tral estimation and so it is not clear if the DM correction 
should be performed. Table [3] indicates these pulsars with a 
'y' or 'n', with the former indicating that we believe that the 
DM correction is likely to improve the residuals. However, 
the DM correction fails to 'improve' PSR J1643— 1224 under 
any metric, even though we measure considerable DM vari- 
ations (see Figure [SJ. As discussed in Section [S] we believe 
that this is due to variations in scattering delay entering the 
DM correction and adding considerable excess noise to the 
corrected residuals. 



6 SCATTERING AND DM CORRECTION 

The most important effect of the ISM on pulsar timing is 
the group delay caused by the dispersive plasma along the 
line of sight. However small scale fiuctuations in the ISM 
also cause angular scattering by a diffractive process. This 
scattering causes a time delay to ~ 0.5^o^/c, where is 
the RMS of the scattering angle and L is the distance to 
the pulsar. This can be significant, particularly at longer 
wavelengths, because it varies much faster with A than does 
the dispersive delay - approximately as A*. In homogeneous 
turbulence one would expect this parameter to be relatively 
constant with time. If so, the delay can be absorbed into the 
pulsar profile and it will have little effect on pulsar timing. 
However if the turbulence is inhomogeneous the scattering 
delay may vary with time and could become a significant 
noise source for pulsar timing. We can study this effect us- 
ing the PPTA pulsar PSR J1939-f 2134. Although this pul- 
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sar is unusual in some respects, the scattering is a property 
of the ISM, not the pulsar, and the ISM in the direction 
of PSR J1939+2134 can be assumed to be typical of the 
ISM in general. PSR J1939+2134 is a very strong source 
and the observing parameters used for the PPTA are well- 
suited to studying its interstellar scattering. The time delay, 
to, can be estimated from the bandwidth of the diffractive 
scintillations, vo, in a dynamic spectrum using the relation- 
ship to = 1/2tivo. In fact it is extremely variable, as can be 
seen in Figure [7] The RMS of to (52 ns at 20 cm) is about 
28% of the mean. We can expect this to increase by a fac- 
tor of (1400 MHz/700 MHz)" = 16 at 50cm cm. Thus in the 
50cm-cm ToAs there will be delays with RMS variations 
of ~ 830 ns, which do not fit the dispersive behaviour. 
This will appear in the estimate of toivi at 20 cm, attenuated 
by a factor of ((1400 MHz/700 MHz) ^ - 1) = 3 (Equation 
[2}. Therefore the DM correction will bring scattering noise 
from the 50-cm band to the 20-cm band with RMS variation 
~ 270 ns, 5.3 times larger than the scattering noise intrinsic 
to the 20-cm observations. This analysis is corroborated by 
the structure function of DM for this pulsar shown in Figure 
[6l which shows a flattening to about 1 /xs^ at small time lags. 
This implies a white process with RMS variations of about 
500 ns, consistent with that expected from scattering. 

We have correlated the variations in ta with the 20-cm 
residuals before correction and find 18% positive correla- 
tion. This is consistent with the presence of 52 ns of com- 
pletely correlated noise due to to added to the ToA measure- 
ment uncertainty of the order of 200 ns. PSR J1939-I-2134 
is known to show ToA vari ations that are correl ated with 
the intensity scintillations (|Cognard et al.1 Il995h but are 
much stronger th an expected for homogeneous turbulence 
l|Coles et aLlboiOl ). Thus we are confident that the observed 
variation in to is showing up in the 20-cm residuals. We 
expect that contribution to increase in the DM corrected 
residuals to about 300 ns. However this is very difficult to 
measure directly because the DM correction is smoothed 
and the 50cm-cm observations are not simultaneous with 
the 20-cm observations. 

This effect increases rapidly with longer wavelength. If 
we had used 80-cm observations for DM correction, the RMS 
at 80 cm would have been ~ 10 /iS and this would have 
been reduced by a factor of 12 to an RMS of 800 ns in the 
corrected 20 cm residua ls. Clearly use of low freque ncy an- 
tennas such as GMRT dJoshi fc Ramakrishnall20(M ) or LO- 
FAR (jStappers et"alll201ll l for correcting DM fluctuations 
in PTAs will have to be limited to weakly scattered pulsars. 
This is an important consideration, but it should be noted 
that the four PPTA pulsars that provide the best timing 
are all scattered much less than J1939-f2134 - all could be 
DM corrected with 80-cm observations or even with longer 
wavelengths. On the other hand there are four PPTA pul- 
sars that are scattered 20 to 80 times more strongly than 
J1939-I-2134 and even correction with 50cm-cm data causes 
serious increases in the white noise. 

An extreme example is PSR J1643— 1224. Under the 
above assumption, the expected white noise (2.0 /xs) due 
to scattering at 20 cm exceeds the radiometer noise (0.63 
^s). The white scattering noise at 50cm cm is much larger 
(32 jis) and about a third of this makes its way into the 
DM-corrected residuals at 20 cm. This is also corroborated 
by the structure function for this pulsar in Figure (6] which 
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Figure 7. Diffractive scattering delay, io, measured from scin- 
tillation bandwidth, vq, in observations of PSR J1939+2134 at a 
wavelength of 20 cm. The error bars are derived from the fit for 
and so are roughly proportional to io. 

shows a flattening to about 10 ^s^ at small lags. This im- 
plies a white process with RMS variation of ~ 3 /is which is 
consistent with that expected from scattering. Indeed, this 
pulsar is the only pulsar with significant DM variations for 
which the DM correction increases the noise in the 20-cm 
residuals under all metrics. It is also important to note that 
observing this source at the same frequencies with a more 
sensitive telescope will not improve the signal to noise ra- 
tio, because the noise, both before and after DM-correction, 
is dominated by scattering. However using a more sensitive 
telescope could improve matters by putting more weight on 
observations at 10 cm, where scattering is negligible. 

Finally however, we note that the usefulness of long 
wavelength observations would be greatly improved if one 
could measure and correct for the variation in scattering de- 
lays. This may be possible using a technique such as cyclic 
spectroscopy, however this has only been done in ideal cir- 
cumstances and with sign al-to-noise ratio such that individ- 
ual pulses are detectable (|Demorestll201lf ). It is still unclear 
if such techniques can be generalised to other observations, 
or if this can be used to accurately determine the unscat- 
tered ToA. 



7 SCHEDULING FOR DM CORRECTION 

If there were no DM variation, one would spend all the ob- 
serving time at the wavelen gth for which the pulsar has the 
greatest ToA precision f see [Manchester et al.|[2"012l for dis- 
cussion on choice of observing wavelength). The reality is, 
of course, that we need to spend some of the time observing 
at other wavelengths to correct for DM variations. In this 
section we will present a strategy for choosing the observing 
time at each wavelength, attempting to optimise the signal 
to noise ratio of the common- mode signal, tcM . We take the 
PPTA observations at Parkes as our example, but this work 
can easily be generalised to any telescope. 

At Parkes it is possible to observe at wavelengths of 10 
and 50cm cm simultaneously because the ratio of the wave- 
lengths is so high that the shorter wavelength feed can be 
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located co-axially inside the longer wavelength feed. How- 
ever the 20-cm receiver does not overlap with either 10 or 
50cm cm and so must be operated separately. 

As noted earlier we can write tcM = Yli ^t^oi so the vari- 
ance of tcM is given by cttcm = Xli bfcri- However, the so- 
lution is not linear in the o"^ terms because they also appear 
in the bi coefficients. At present, observing time at Parkes is 
roughly equal between the two receivers, so we use the exist- 
ing power spectral densities as the reference. We will assume 
that the total observing time is unity and the time devoted 
to 10 and 50cm cm, which are observed simultaneously, is 
< a < 1. 

The variance of any toi is inversely proportional to the 
corresponding observing time. We use to, 20 as the refer- 
ence because it usually has the smallest TOA uncertainty 
in PPTA observations. Therefore, we define (j|q = 1/(1 — a) 
as the reference. Then we assume that, with equal observing 
time (Tio = xa2o and aso = 2/0-20, so as scheduled we would 
have (Tio = x/a and (T50 = y/a. We can then determine the 
increase in white noise caused by correcting for dispersion as 
a function of a, the time devoted to 10 and 50cm cm. The re- 
sults are shown for all the PPTA pulsars in Table [T] One can 
see that all the pulsars are different and the optimal strate- 
gies range from a « 0.2 to a « 1.0 (i.e. 100% of time spent 
using the dual 10 and 50cm cm system). For the four 'best' 
pulsars, PSRs J0437-4715, J1713-h0747, J1744-1134, and 
J1909-3744, the optimal strategy has a > 0.7. 

This suggests that a useful improvement to PTA perfor- 
mance could come from deploying broadband receivers, so 
that correction for DM(t) can be done with a single observa- 
tion. This also has the benefit of reducing the difficulties of 
aligning pulsar profiles measured with different receivers at 
different times, and would therefore allow for more accurate 
measurement of DM variations. 



8 THE INTERSTELLAR MEDIUM 

The PPTA DM(t) observations provide an interesting pic- 
ture of the ionised ISM on au spatial scales. The overall 
picture can be seen in Figures [5] and [6] In Figure [5] it is ap- 
parent that 17 of the 20 PPTA pulsars have measurable 
DM(t) variations. In Figure [H] it can be seen that 13 of 
these 17 show power-law structure functions, as expected in 
the ensemble average Of these, eight are roughly consistent 
with an extrapolation from the diffractive scales at the Kol- 
mogorov spectral exponent, an average dynamic range of 4.8 
decades. However five are considerably higher than is pre- 
dicted by a Kolmogorov extrapolation. They may be locally 
Kolmogorov, i.e. an inner scale may occur somewhere be- 
tween the diffractive scale and the directly measured scales 
of 100 to 2000 days, but establishing this would require a 
detailed analysis of the apparent scintillation velocity which 
is beyond the scope of this paper. Two of these five pulsars, 
J1045-4509 and J1909-3744, were already know n to be in- 
consi stent with a Kolmogorov spectral exponent (|You et al.l 
120071 ). and it is clear, with the additional data that are now 
available, that J1024-0719, J1643-1224 and J1730-2304 
should be added to this list. When the spatial power spec- 
trum of a stochastic process is steeper than the Kolmogorov 
power-law, it can be expected to be dominated by linear 
gradients and show an almost quadratic structure function. 



Indeed inspection of Figure [5] shows that the 5 steep spec- 
trum pulsars all show a strong linear gradient in DM(t). 

The time series DM(t) shown in Figure [5] often show 
behaviour that does not look like a homogeneous stochastic 
process. For example, PSR J1603— 7202 shows a large in- 
crease for ~250 days around MJD 54000 and J0613-0200 
shows clear annual modulation. The increase in DM for 
J1603— 7202 suggests that a blob of plasma moved through 
the line of sight. If we assume the blob is halfway between 
the pulsar and the Earth, the line of sight would have moved 
by about 0.5 au in this time, and if the blob were spheri- 
cal it would need a density of ~200 cm~^. This value is 
high, but comparable to other density estimat es for au-scale 
structure based on 'extrem e scattering events' (|Fiedler et al.l 
1 19871 : ICognard et allll993l ) . 

We computed the power spectra of DM(t) for all the 
pulsars to see if the annual modulation that is clear by 
eye in PSR J0613— 0200 is present in any of the other pul- 
sars. For four pulsars we find a significant (> 5-a) detection 
of an annual periodicity, PSRs J0613-0200, J1045-4509, 
J1643-1224 and J1939-f2134. 

The most likely explanation for the annual variation in 
DM(t) is the annual shift in the line of sight to the pulsar 
resulting from the orbital motion of the Earth. The trajec- 
tory of the line of sight to three example PPTA pulsars are 
shown in Figure [S] The relatively low proper motion and 
large parallax of the PPTA pulsars means that the trajec- 
tory of the line of sight to many of the PPTA pulsars show 
pronounced ripples. However, unless the trajectory is a tight 
spiral, the annual modulation will only be significant if there 
is a persistent gradient in the diffractive phase screen. 

The presence of persistent phase gradients and annual 
modulation in J1045-4509 and J1643-1224 is not surpris- 
ing because the ISM associated with each of these pulsars 
has a steeper than Kolmogorov power spectrum. Indeed, 
the measured DM(t) for these pulsars do show a very lin- 
ear trend, which in itself evidence for a persistent phase 
gradient. The other steep spectrum pulsars, J1024— 0719, 
J1730-2304 and J1909-3744, have higher proper motion, 
which reduces the amplitude of the annual modulation rela- 
tive to the long term trend in DM(t). We note that the spec- 
tral analyses for PSRs J1024-0719 and J1909-3744 suggest 
annual periodicities, and it may be possible to make a sig- 
nificant detection by combining the PPTA data with other 
data sets. 

PSR J1939-I-2134 does not show a steep spectrum, how- 
ever its proper motion is very low compared to its parallax, 
and therefore the trajectory spirals through the ISM, reduc- 
ing the requirement for a smooth phase screen. The annual 
modulation of J0613— 0200 may be somewhat different, since 
it does not have a steep spectrum and although the proper 
motion is small the trajectory does not spiral (see Figure 
[8|. This suggests that for J0613-0200 the turbulence could 
be anisotropic with the slope of the gradient aligned with 
the direction of the proper motion. Anisotro pic structures 
are believed to be quite c ommon in the ISM jCordes et al.l 
I2OO6I : iBrisken et alTl2010l V However one can imagine vari- 
ous other ways in which this could occur, particularly in 
an inhomogeneous random process, and inhomogeneous tur- 
bulence on an a u spatial scale is also believed to be com- 
mon in the ISM ([Stinebring et al.ll200ll : ICordes et al 1 l2006l : 
iBrisken et alll201of r 
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Persistent spatial gradients will cause a refractive shift 
in the apparent position of the pulsar, and because of dis- 
persion the refraction angle will be wavelength dependent. 
This refractive shift appears in the timing residuals as an 
annual sine wave which changes in amplitude like A'^. When 
the DM(t) is corrected this sine wave disappears and the in- 
ferred position becomes the same at all wavelengths. These 
position shifts are of order 10"" (A/20 cm) arcseconds for 
all four pulsars. 

Note that the trajectory of the lines of sight shown on 
Figure [8] may appear quite non-sinusoidal, but the annual 
modulation caused by the Earth's orbital motion in a linear 
phase gradient will be exactly a sine wave superimposed on a 
linear slope due to proper motion. This will not generate any 
higher harmonics unless the structure shows significant non- 
linearity on an au scale. We do not see second harmonics of 
the annual period, which suggests that the spatial structure 
must be quite linear on an au scale. 

Annual variations in DM are also observed in pulsars 
for which the line of sight passes close to the Sun because of 
free e lectrons in the solar wind (|Ord et al.ll2007l : IYou et"al] 
I2OI2I ). In the PPTA, a simple symmetric model of the solar 
wind is used to remove this effect, but this is negligible for 
most pulsars. For the three pulsars where it is not negligible, 
the effect of the solar wind persists only for a few days at the 
time when the line of sight passes closest to the Sun. Neither 
the magnitude, phase nor shape of the variations seen in our 
sample can be explained by an error in the model of the 
solar wind. Changes in ionospheric free electron content can 
be ruled out for similar reasons. 

In summary the ISM observations are, roughly speak- 
ing, consistent with our present understanding of the ISM. 
However the data will clearly support a more detailed anal- 
ysis, including spectral modelling over a time scale range 
in excess of 10^ from the diffractive scale to the duration 
of the observations. It may also be possible to make a 2- 
dimensional spatial model of the electron density variations 
for some of the 20 PPTA pulsars. Although this would be 
useful for studying the ISM and in improving the DM cor- 
rection, such detailed modelling is beyond the scope of this 
paper. 



9 CONCLUSIONS 

We find that it is necessary to approach the problem of es- 
timating and correcting for DM(f) variations iteratively, be- 
ginning with a pilot analysis for each pulsar and refining that 
analysis as the properties of that pulsar and the associated 
ISM become clearer. Each pulsar is different and the ISM in 
the line of sight to each pulsar is different. The optimal anal- 
ysis must be tailored to the conditions appropriate for each 
pulsar and according to the application under consideration. 

We sample the DM(t) just often enough that the varia- 
tions in DM are captured with the minimum amount of addi- 
tional white noise. Likewise, we must also sample a common- 
mode signal tcM{t) at the appropriate rate. In this way we 
can correct for the DM variations at frequencies where it is 
necessary, and we can include tcM{t) at frequencies where 
it is necessary, but not fit for either at frequencies where the 
signal is dominated by white noise. 

By including the common-mode signal in the analysis 
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Figure 8. Trajectories through the ISM of the line of sight to 
PSRs J0613-0200 (dashed Hne), J1643-1224 (soUd black line) 
and J1909— 3744 (grey line). It was assumed that the scattering 
takes place half way between the pulsar and the Earth and the 
motion of the plasma was neglected. The trajectories are marked 
with a cross at the DM sampling interval of 0.25 yr. 



we preserve the wavelength-independent signals of interest 
for pulsar timing arrays and we improve the estimate of 
the pulsar period and period derivative. Without estimat- 
ing the common mode, a significant fraction of wavelength- 
independent signals, such as: errors in the terrestrial clocks; 
errors in the planetary ephemeris; and the effects of gravita- 
tional waves from cosmic sources, would have been absorbed 
into the DM correction and lost. 

We have applied this technique to the PPTA data set, 
which improves its sensitivity for the detection of low fre- 
quency signals significantly. The estimated DM{t) also pro- 
vides an unparallelled measure of the au scale structure of 
the interstellar plasma. In particular it confirms earlier sug- 
gestions that the fiuctuations often have a steeper than Kol- 
mogorov spectrum, which implies that an improved physical 
understanding of the turbulence will be necessary. We also 
find that persistent phase gradients over au scales are rel- 
atively common and are large enough to cause significant 
errors in the apparent positions of pulsars unless DM cor- 
rections are applied. 
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APPENDIX A: CONSTRAINED LEAST 
SQUARES FITTING IN TEMP02 

The least squares problem of fitting the timing model to the 
residuals can be written in matrix form as 

R = MP + E. (Al) 



Here R is a column vector of the timing residuals, P is 
a column vector of fit parameters, including DM(tj) and 
icM(ij) as well as the other timing model parameters. M 
is a matrix describing the timing model and E is a column 
vector of errors. The least-squares algorithm solves for P, 
matching MP to R with a typical accuracy of E. 

The sampled time series DM(tj) and icM(ij) are co- 
variant with the timing model, so they must be con- 
strained to eliminate that covariance or the least squares 
solution will fail to converge on a unique solution. These 
constraints have the form of linear equations of DM(tj) 
and tcmitj), such as: Y^DM{tj) = 0; Y^tcM{tj) = 0; 
EijicMfe) = 0, EijicMfe) = 0; EsinMi)tcMfe) = 0; 
'^cos{ujtj)tcyi{tj) = 0; etc. Augmented with these equa- 
tions, the least-squares problem becomes 



R 




M 


P + 


E 


C 




B 


e 



where B is a matrix describing the constraints, e is a column 
vector of weights for the constraints. In our case C = 0, 
though it need not be in general. The least-squares solution 
will then find a vector P that matches both MP to R, with 
a typical accuracy of E, and also matches BP to C, with a 
typical accuracy of e. By making e very small we can enforce 
the constraints with high accu racy. This scheme has b een 
called 'the method of weights' (jColub fc van Loanlll996l ). 

If the uncertainties in the estimates of DM(tj) and 
tcM{tj) are not expected to be equal, for instance if the 
different observing wavelengths are irregularly sampled and 
the To A uncertainties are variable across sampling windows, 
then it can be advantageous to use weighted constraints. 
Then the constraints take the form ^ WjDM{tj) = 0, and 
we need to estimate the uncertainties of the parameters to 
obtain the optimal weights. These uncertainties can be de- 
termined from the least-squares solution in which the timing 
residuals are described purely by Equation Q . This problem 
is linear and the covariance matrix of the parameters can be 
written in closed form without even solving for the parame- 
ters. The diagonal elements of the covariance matrix are the 
variances of the parameters and the weights, Wj, are the 
inverse of the square roots of the corresponding variances. 



APPENDIX B: RELATION BETWEEN THE 
STRUCTURE FUNCTION AND 
POWER-SPECTRAL DENSITY 

The structure function D{t), of a time series y{t), is well 
defined if y{t) has stationary differences 

D{r) = {{y{t)^y{t + r)f). (Bl) 

If y{t) is wide-sense stationary D{t) can be written in terms 
of the auto covariance C(r) by expansion of the square 

D{r) = {y(tf) + {y(t + rf)-2{y{t)y(t + r)) 

= 2(C(0)-C(r)) (B2) 

If y{t) is real valued then by the Wiener-Khinchin theorem, 

C{r)= rcos{2nfr)P{f)df, (B3) 
Jo 

where P(f) is the one-sided power spectral density of y(t). 
Thus we can then write the structure function in terms of 
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the power spectral density as 

D{t) = A 2(l-cos(2^/r))P(/)d/, (B4) 
Jo 

It should be noted that this expression for D{t) is valid if 
D{t) exists. It is not necessary that C(r) exist. For the case 
of a power-law, P{f) — Af^", we can change variables using 
X = fr, and obtain 

D(t)=t''~^A 2(1 -cos(27r2;))2;"''da:. (B5) 
Jo 

The integral (Int) above converges if 1 < q < 3, yielding 

Int. = 2"7r""^ sin(-a7r/2)r(l - a), (B6) 

where F is the Gamma function. Thus for Kolmogorov tur- 
bulence, with exponent a — 8/3, we have Int ~ 89.344 and 
the power spectrum can be written 

P{f) ~ 0.0112D(r)r~^''V"*'^^- (B7) 
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ABSTRACT 

Signals from radio pulsars show a wavelength-dependent delay due to dispersion in 
the interstellar plasma. At a typical observing wavelength, this delay can vary by tens 
of microseconds on five-year time scales, far in excess of signals of interest to pulsar 
timing arrays, such as that induced by a gravitational- wave background. Measurement 
of these delay variations is not only crucial for the detection of such signals, but 
also provides an unparallelled measurement of the turbulent interstellar plasma at au 
scales. 

In this paper we demonstrate that without consideration of wavelength- 
independent red-noise, 'simple' algorithms to correct for interstellar dispersion can 
attenuate signals of interest to pulsar timing arrays. We present a robust method for 
this correction, which we validate through simulations, and apply it to observations 
from the Parkes Pulsar Timing Array. Correction for dispersion variations comes at 
a cost of increased band-limited white noise. Wc discuss scheduling to minimise this 
additional noise, and factors, such as scintillation, that can exacerbate the problem. 

Comparison with scintillation measurements confirms previous results that the 
spectral exponent of electron density variations in the interstellar medium often ap- 
pears steeper than expected. We also find a discrete change in dispersion measure of 
PSR J1603-7202 of ~ 2 x lO'^ cm'^pc for about 250 days. We speculate that this has 
a similar origin to the 'extreme scattering events' seen in other sources. In addition, 
we find that four pulsars show a wavelength-dependent annual variation, indicating 
a persistent gradient of electron density on an au spatial scale, which has not been 
reported previously. 

Key vifords: pulsars: general — ISM: structure — methods: data analysis 



1 INTRODUCTION 

The fundamental datum of a pulsar timing experiment is 
the time of arrival (ToA) of a pulse at an observatory. In 
practise, the ToA is referred to the solar-system barycen- 

* Email: mkeith@pulsarastronomy.net 



tre in a standard time frame (e.g., barycentric coordinate 
time). This barycentric arrival time can be predicted using 
a 'timing model' for the pulsar. The difference between the 
barycentric ToAs and the arrival times predicted by the tim- 
ing model are termed residuals. The timing model can be re- 
fined using a least-squares fitting procedure to minimise the 
residuals, as performed by, e.g., the Tempo2 software (?). 
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Since the timing model is always incomplete at some level, 
we always see some level of post-fit residuals, which are typi- 
cally a combination of 'white' noise due to the uncertainty in 
the ToA measurement and 'red' (i.e., time-correlated) sig- 
nal. For the majority of known pulsars the dominant red 
signal is caused by the intrinsic instability of the pulsar, 
and termed 'timing noise' (e.g., ?). However, the subset of 
millisecond pulsars are stable enough that other red signals 
are potentially measurable (?). Pulsar timing array projects, 
such as the Parkes Pulsar Timing Array (PPTA; ?), aim to 
use millisecond pulsars to detect red signals such as: errors 
in the atomic time standard (?); errors in the Solar System 
ephemeris (?); or the effect of gravitational waves (???). 
Each of these signals can be distinguished by the spatial 
correlation, i.e., how pulsars in different directions on the 
sky are affected. However, at typical observing wavelengths 
and time-spans, the variation of the dispersive delay due to 
turbulence in the ionised interstellar medium (ISM) domi- 
nates such signals (?). Fortunately for pulsar timing exper- 
iments, these delays can be measured and corrected using 
observations at multiple wavelengths. 

The dispersive group delay is given by 



t-DM = A 



2-KmeC?' 



ne{l)dl 

path 



(1) 



where A is the barycentric radio wavelengtfQ. The path in- 
tegral of electron density is the time-variable quantity, fn 
pulsar experiments this is termed 'dispersion measure', DM, 
and given in units of cm~''pc. In principle, the instantaneous 
DM can be computed from the difference of two arrival times 
from simultaneous observations at different wavelengths, or 
more generally by fitting to any number of observations at 
more than one wavelength. 

The question of estimation and correction of DM{t) has 
previously been considered by ?. They chose a 'best' pair 
of wavelengths from those available and estimated the DM 
at every group of observations. These observation groups 
were selected by hand, as was the choice of wavelengths. 
Regardless of how the analysis is done, the estimated DM 
always contains white noise from differencing two observa- 
tions, and correcting the group delay always adds that white 
noise to the arrival times. However the DM{t) variations are 
red, so they only need to be corrected at frequencies be- 
low the 'corner frequency' at which the power spectrum of 
the DM-caused fluctuations in group delay is equal to the 
power spectrum of the white noise in the T)M{t) estimate. 
To minimise the additional white noise, they smoothed the 
DM(t) estimates over a time Ts to create a low-pass filter 
which cuts off the DM variations, and the associated white 
noise, at frequencies above the corner frequency. In this way, 
they avoided adding white noise at high frequencies where 
the DM-correction was unnecessary. Of course the added 
'white' noise is no longer white; it is white below the corner 
frequency, but zero above the corner frequency. 

Here we update this algorithm in two ways. We use 
all the observed wavelengths to estimate DM(t) and we in- 
tegrate the smoothing into the estimation algorithm auto- 



^ To avoid confusion, in this paper we will use wavelength for 
the radio wavelength and frequency to describe the fluctuation of 
time variable processes. 



matically. Thus, the algorithm can easily be put in a data 
'pipeline'. We show the results of applying this new algo- 
rithm to the PPTA data set, which is now about twice as 
long as when it was analysed by ?. Additionally, we demon- 
strate that our algorithm is unbiased in the presence of 
wavelength-independent red signals, e.g., from timing noise, 
clock error, or gravitational waves; and we show that failure 
to include wavelength-independent red signals in the esti- 
mation algorithm will significantly reduce their estimated 
amplitude. 



2 THEORY OF DISPERSION REMOVAL 

We assume that an observed timing residual is given by 
ioBS = ^CM -I- tDM(A/AREF)^ where tcm is the common- 
mode, i.e., wavelength-independent delay and tuM is the dis- 
persive delay at some reference wavelength Aref. Then with 
observations at two wavelengths we can solve for both tcM 
and toM- 



ioM = (^OBS.l — iOBS,2)AREF/(Al 



icM = (tOBS,2Ai — tOBS,lA2)/(Al 



Ai), 



Ai). 



(2) 
(3) 



In a pulsar timing array, tcM would represent a signal of 
interest, such as a clock error, an ephemeris error, or the ef- 
fect of a gravitational wave. The dispersive component tDM 
would be of interest as a measure of the turbulence in the 
ISM, but is a noise component for other purposes. It is im- 
portant to note that tuM is independent of tcM so one can 
estimate and correct for the effects of dispersion regardless 
of any common- mode signal present. In particular, common- 
mode red signals do not cause any error in toM- 

If more than two wavelengths are observed, solving for 
tcM and tuM becomes a weighted least-squares problem, and 
the standard deviation of the independent white noise on 
each observation is needed to determine the weighting fac- 
tors. For wavelength i, we will denote the white noise by 
tw,i and its standard deviation by at so the observed timing 
residual is modelled as 



tOBS.i = tCM + toui^i / ^REf)^ + tw,i 



(4) 



The weighted least-squares solutions, which are minimum 
variance unbiased estimators, are 



toM = Aref ^ ^ ioBS,iAi /(Tj^ — 

i i 

^A?/a?^toBS../a?)/A (5) 

i i 
i i 

^A?/a?^toBS,.A?/a?)/A. (6) 

i i 

Here A is the determinant of the system of equations, 

i i i 

If one were to model only the dispersive term toM, the 
weighted least-squares solution would become 

I]i*OBS,iAi/crf 



toM = Aj, 



E.At/a? 



(7) 
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However if a common- mode signal is present, this solution 
is biased. The expected value is 

(ioivi) ~ toM + ^cmAref /J, 2 ■ (8) 

Some of the 'signal' tcm is absorbed into toM reducing the 
effective signal-to-noise ratio and degrading the estimate of 
DM. We will demonstrate this bias using simulations in Sec- 
tion [1 

It is important to note that the dispersion estimation 
and correction process is linear - the estimators toM and 
tcM are linear combinations of the residuals. The corrected 
residuals toss, cor, i = ioBS,i — (Ai/AREF)^toM, are also lin- 
ear combinations of the residuals. We can easily compute 
the white noise in any of these quantities from the white 
noise in the residuals. For example, we can collect terms 
in Equations ([5]) and ([6]) obtaining turn = X^i ciitoBS.i and 
icM = J2i bitoBS,i, where 

a, = AW(A?/<T?^l/a|-l/a?^A,Va|)/A (9) 

= (l/af^A|/a|-A?/a?^A,Va|)/A. (10) 

j j 

Then, the white noise variances of the estimators can be 
written as a^uM = Ei and cttcm = J2i blcrf. 

The actual PPTA observations are not simultaneous at 
all frequencies, so we cannot normally apply Equations (O 
and © directly (?). We discuss how the least squares solu- 
tions for toM and tcM can be obtained by including them in 
the timing model in the next section. However it is useful to 
have an analytical estimate of the power spectral density of 
the white noise that one can expect in these estimators and 
in the corrected residuals. At each wavelength \i we have a 
series of Ni error estimates cr^j . The variance of the weighted 
mean is aj^i = 1/ Ej ^/'^ij- This is the same as if we had a 
different number A'^ of observations at this wavelength each 
of variance = a^^N. Thus, for planning purposes we can 
compute ami for each wavelength and conceptually resample 
each wavelength with an arbitrary number (A) of samples. 
Equations ((5]), ((6]), ((9]), and (fTO)) are invariant under scaling 
of all o"i by the same factor so one can obtain the coefficients 
Ui and bi using ami in place of ai so the actual number (Ai) 
of samples need not enter the equations. 

If one had a series of N samples over a time span of 
ToBS each with variance a^, the spectral density of the white 
noise would be = 2roBS a^ /N = 2T'obs am- We can ex- 
tend this to a weighted white noise spectral density using 
the variance of the weighted mean. So the power spectral 
densities Piu,i play the same role as a^ in Equations ((S]), 
((G]), (O and (|10p . The coefficients {ai} and {bi} are func- 
tions of Ai and Pw,i - Then we find P,„,tdm ~ Ei "^Pwa and 
P-w,Tcm = Ei biPw,i. 

Perhaps the most important property of these estima- 
tors is that Pw,TCM is less than or equal to the white noise 
spectrum of the corrected residuals Pw,cor,i in any band. 
Equality occurs when there are only two wavelengths. The 
values of Pm,i, Pw,cor,i, Pm,TDM and P™,TCM are given for 
the PPTA pulsars in Table [T] Here Pw,tdm is given at the 
reference wavelength of 20 cm. 

The situation is further complicated by red noise which 
depends on wavelength, but not as A^. For example, diffrac- 



tive angular scattering causes variations in the group delay, 
which scale as the scattered pulse width, i.e. approximately 
as A" (?). Clearly such noise will enter the DM correction 
process. It can have the unfortunate effect that scattering 
variations, which are stronger at long wavelengths, enter the 
short wavelength corrected residuals even though they are 
negligible in the original short wavelength data. This will be 
discussed in more detail in Section [G] 



3 DISPERSION CORRECTION TECHNIQUE 

Rather than solving for tcM and toM for every group of 
observations, or re-sampling observations at each wave- 
length to a common rate, it is more practical to include 
parametrised functions for tcM{t) and DM{t) in the tim- 
ing model used to obtain the timing residuals. To provide a 
simple and direct parametrisation we use piece-wise linear 
models defined by fixed samples tcM{tj) and DM(tj) for j = 
1,...,A.. 

It is also required to introduce some constraints into the 
least-squares fitting to prevent covariance with other model 
parameters. For example, the values of DM{tj) axe natu- 
rally covariant with the mean dispersion measure parame- 
ter, DMo, which is central to the timing model. To eliminate 
this covariance, we implement the linear equality constraint 
that J2i=i'DM{tj) = 0. Additionally, the series tcM(ij) is 
covariant with the entire timing model, however in practise 
the sampling interval is such that it responds very little to 
any orbital parameters (in the case of binary systems). We 
constrain tcuitj) to have no response to a quadratic poly- 
nomial, or to position, proper motion, and parallax. These 
constraints are implemented as part of the least-squares fit 
in Tempo2, as described in Appendix 1X1 

The choice of sampling interval, Tg is essentially the 
same as in ?. The process of fitting to a piece- wise linear 
function is equivalent to smoothing the T)M{t) time series 
with a triangle function of base 2Ta. This is a low pass filter 
with transfer function Htn{f) = {sm{nfTs)/TrfTs)^. We ad- 
just Ts such that the pass band approximately corresponds 
to the corner frequency fc at which the power spectrum 
of the DM delays, Ptdm, exceeds that of the white noise, 
Pw, TDM- Note that this corner frequency is independent of 
reference wavelength at which turn is defined. 

To determine this corner frequency we need an estimate 
of the power spectrum of tnu, so the process is inherently 
iterative. We can obtain a first estimate of Ptdm(/) from the 
diffractive time scale, Tdis, at the reference wavelength. For 
signals in the regime of strong scattering, which includes all 
PPTA observations, rdifi is the time scale of the diffractive 
intensity scintillations. For the PPTA pulsars, Tdig is usually 
of the order of minutes and can be estimated from a dynamic 
spectrum taken during normal observations (see e.g. ?). 

Rather than directly compute Ptdm, it is attractive to 
begin with the structure function, which is a more conve- 
nient statistic for turbulent scattering processes and is more 
stable when only a short duration is available. The structure 
function of tuM is given by 

I>TDM(r) = {{tuM{t) - tuM{t + T)f} = (\/2ncf D^{t) , 

(11) 

where D^ij) is the phase structure function. If we assume 
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Table 1. The estimated power spectral density before {Pw) and after (P^j.cor) correction of the white noise for each PPTA pulsar at 
each of the three wavelengths, and the expected white noise power spectral density in the 'common mode' signal (P^^tcm) ^'^d in (dm 
at 20 cm (Pii,,tdm)i all expressed relative to the power spectral density of the uncorrected 20-cm residuals. Also shown is the effect of 
optimising the observing time, expressed as the ratio of Pu.'.TCM estimated for optimal observing and P^j^tcm with the cmrent observing 
strategy (q = 0.5), and Qopt the optimal fraction of time spent using the dual 10- and 50cm-cm observing system. 
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0.6 
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33 


12 


33 


2.9 


3.2 


2.9 


1.4 


1.0 


0.5 


J1045-4509 
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17 


3.1 


17 


2 


2 


1.9 


0.43 


0.9 


0.3 


J1600-3053 
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15 


3.3 


2.1 


4.6 


1.8 


1.2 


0.9 
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1.5 
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0.3 
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1.0 
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0.94 


0.84 


0.22 


0.8 
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10 
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1.2 


1.0 


0.5 


J1744-1134 
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12 


4.4 


2.2 


3.7 


1.9 


1.1 


0.9 


0.7 


J1824-2452A 
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30 


40 


32 


5.5 


8.7 


5.4 


4.1 


0.9 


0.7 


J1857+0943 
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20 
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0.9 


0.6 
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4.3x10-30 


0.51 
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1.1 
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0.5 
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0.6 


1.0 
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14 
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2.2 


2.1 
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1.0 


0.4 
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19 


2 


1.9 


1.9 


0.4 


0.9 


0.3 
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2.5 
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2.1 


2.1 


2.1 


0.39 


0.8 


0.3 
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8.2 


15 


8.8 


2.7 


4 


2.6 


1.5 


1.0 


0.6 



that the electron density power spectrum has an expo- 
nent of -11/3, i.e., Kolmogorov turbulence, then D^{t) = 
(r/rdiff)^^^ (?). The structure function Dtdm(t") can there- 
fore be estimated from rdift, or directly from the toM{t) once 
known. 

As described in Appendix |B] we can use the structure 
function at any time lag r to obtain a model power spectrum 
using 

Ptdm(/) ^ 0.0112 DTDM(r)r-'''='(spy)-i/V"*^' (12) 

The term (spy) is the number of seconds per year. Here 
Dtom is in s2, r in s, / is in yr"""^ and Ptdm is in yr^. 

The spectrum of the white noise can be estimated from 
the ToA measurement uncertainties as discussed in section 
2. However, often there are contributions to the white noise 
that are not reflected in the measurement uncertainties and 
so we prefer to estimate directly from the power spec- 
trum of the residuals. 



4 TEST ON SIMULATED OBSERVATIONS 

When dealing with real data sets it is not trivial to show that 
the DM-corrected residuals are 'improved' over simply tak- 
ing residuals from the best wavelength (?). This is because 
much of the variations in DM are absorbed into the fit for 
the pulsar period and period derivative. Therefore the root- 
mean-square deviation (RMS) of the residuals from a single 
wavelength may not decrease significantly even though the 
RMS of the DM(f) variations that were removed is large. To 
demonstrate that the proposed procedure can estimate and 
remove the dispersion, and that it is necessary to include 



the common-mode in the process, we perform two sets of 
simulations. 

The observing parameters, i.e., Tobs, Ni, aij, -Ddm(t), 
of both simulations are based on the observations of PSR 
J1909-3744 in the PPTA 'DRl' data set (?). We find it 
useful to demonstrate the performance of the DM correc- 
tion process in the frequency domain, but it is difficult to 
estimate power spectra of red processes if they are irregu- 
larly sampled. Therefore we first use simulations of regularly 
sampled observations with observing parameters similar to 
those of PSR J1909— 3744 to demonstrate the performance 
of the DM correction algorithm. Then we will simulate the 
actual irregularly sampled observations of PSR J1909— 3744 
to show that the ultimate performance of the algorithm is 
the same as in the regularly sampled case. 

4.1 Regular sampling, equal errors 

We will compare the power spectra produced after fitting 
for DM(t) with and without simultaneously fitting for a 
common-mode signal. To generate the simulated data sets, 
we first generate idealised ToAs that have zero residual from 
the given timing model. Then we add zero-mean stochastic 
perturbations to the ideal ToAs to simulate the three compo- 
nents of the model: (1) independent white noise, correspond- 
ing to measurement error; (2) wavelength independent red 
noise, corresponding to the common-mode; (3) wavelength 
dependent red noise representing DM{t). 

We simulate the measurement uncertainty with a white 
Gaussian process, chosen to match the high frequency power 
spectral density of the observed residuals. The simulated 
P™ is 2.2 x 10"3°,4.3 x lO'^o and 2.6 x 10"^^ yr^ at 10, 
20 and 50cm cm respectively. For the common mode we 
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choose a Gaussian process with a spectrum chosen to match 
a common model of the incoherent gravitational wave back- 
ground (GWB), i.e. Pgwb(/) = (AI^^/Utv^)/-''^' (??). 
For the DM we use a Gaussian process with a power spec- 
trum Pdm(/) = ^dm/"*''^, where ^dm is chosen to match 
the observed DM fluctuations in PSR J1909-3744 shown 
in Figures [5] and (6] and the spectral exponent is chosen to 
match that expected for Kolmogorov turbulence (?). The 
levels of Ptdm and Pgwb are similar so that the same sam- 
ple intervals can be used for both T)M{ti) and tcuiti), but 
this is not necessary in general and will not always be desir- 
able. 

For both algorithms, we estimate the pre- and post- 
correction power spectra of the 20-cm residuals in four noise 
regimes: P^, ; P,„ +Pdm; Pw + Pgwb ; and P,„ -f Pdm + Pgwb • 
In order to minimise the statistical estimation error, we aver- 
age together 1000 independent realisations of the spectra for 
each algorithm. We note that although the averaged power 
spectra suggest that the input red noise signals are large, 
the noise on a single realisation is such that the red signals 
are at the limit of detection. To illustrate this, the 90% con- 
fidence limits for both the 1000 spectrum average and for a 
single realisation, are shown on the power spectra in Figures 
[T]and[l 

We show the effect of using the interpolated model for 
DM(t), but not fitting for the common-mode signal tcM(i), 
in Figure [T] This algorithm is well behaved when the GWB 
is not present, as shown in the two lower panels. In this 
case the DM correction algorithm removes the effect of the 
DM variations if they are present and increases the white 
noise below the corner frequency by the expected amount. 
Importantly, when the model GWB is included, i.e., in the 
two top panels, a significant amount of the low-frequency 
GWB spectrum is absorbed into the DM correction. This is 
independent of whether or not DM variations are actually 
present because the DM correction process is linear. 

We show the full algorithm developed for this paper, 
using interpolated models for both T>M{t) and the common- 
mode signal tcM{t), in Figure O One can see that the algo- 
rithm removes the DM if it is present, regardless of whether 
the GWB is present. It does not remove any part of the 
GWB spectrum. When the GWB is not present, as shown 
in the two lower panels, the algorithm remains well behaved. 
As expected, it increases the white noise below the corner 
frequency by a larger factor than in the previous case. This 
is the 'cost' of not absorbing some of the GWB signal into 
the DM correction. Although it has a higher variance than 
for the previous case, our DM{t) is the lowest variance un- 
biased estimator of the DM variations in the presence of 
wavelength-independent red noise. This increase in white 
noise is unavoidable if we are to retain the signal from a 
GWB, or indeed any of the signals of interest in PTAs. 

The power spectra presented in Figure [2] demonstrate 
that the algorithm is working as expected, in particular that 
it does not remove power from any wavelength-independent 
signals present in the data. We note, however, two limita- 
tions in these simulations: the regular sampling and equal 
errors are not typical of observations, nor have we shown 
that the wavelength-independent signal in the post-fit data 
is correlated with the input signal (since our power spec- 
trum technique discards any phase information) . These lim- 
itations will be addressed in the next section. 



4.2 Irregular sampling, Variable error bars 

In order to test the algorithm in the case of realistic sampling 
and error bars, we repeated the simulations using the actual 
sampling and error bars for pulsar J1909— 3744 from the 
PPTA. We use the same simulated spectral levels for the 
GW and DM as in the previous section. The results are also 
an average of 1000 realisations. 

As a direct measure of performance in the estimating 
DM(t), we compute the difference between the DM esti- 
mated from the fit to the residuals, DMcst(t), and the DM 
input in the simulation, DMin(t). To better compare with 
the timing residuals, we convert this error in the DM into 
the error in tumit) at 20cm using Equation ([T]). Note that, 
although the residuals were sampled irregularly, the original 
DMin(t) was sampled uniformly on a much finer grid. Fur- 
thermore, the estimated DMeat (t) is a well defined function 
that can also be sampled uniformly. Thus it is easy to com- 
pute the average power spectrum of this error in tuM{t) as is 
shown in Figure [31 We also plot the spectrum of the initial 
white noise, and the spectrum of the white noise after cor- 
rection. If the algorithm is working correctly the white noise 
after correction should exactly equal the error in toM{t) plus 
the white noise before correction, so we have over plotted the 
sum of these spectra and find that they are identical. 

The spectrum of the error in tDM(i) shows the expected 
behaviour below the corner frequency. Above the corner fre- 
quency (where the correction is zero), it falls exactly like the 
spectrum of t-ouit) itself, i.e., as f~^^^- By comparing the 
right and left panels one can see that the DM correction is 
independent of the GWB. 

We can also demonstrate that the model GWB sig- 
nal is preserved after DM correction, by cross-correlation 
of the input model GWB with the post-correction residu- 
als. If the GWB signal is preserved this cross-correlation 
should equal the auto-correlation of the input GWB signal. 
We show the auto-correlation of the input and four different 
cases of the cross-correlation of the output in Figure [4] The 
cross-correlations are for two bands (20 and 50cm cm shown 
solid and dashed respectively), and for two different fitting 
algorithms (with and without tcM{t) shown heavy and light 
respectively). Again it can be seen that, without fitting for 
the common-mode tcM{t), a significant portion of the GWB 
is lost. In fact, it is apparent from the large negative correla- 
tion at 50cm cm that the 'lost' power is actually transferred 
from the 20-cm residuals to those at 50cm cm. Although it 
may be possible to recover this power post-fit, it is not clear 
how to do this when the GWB and DM signals are unknown. 
Finally, we note that when the common mode is used, the 
50cm-cm residuals preserve the GWB just as well as the 20- 
cm residuals, even though they carry the majority of the 
DM(t) variation. 



4.3 The Robustness of the Estimator 

The proposed DM correction process is only optimal if the 
assumptions made in the analysis are satisfied. The pri- 
mary assumptions are: (1) that there is an unmodelled 
common-mode signal in the data; (2) that the residuals 
can be modelled as a set of samples toi{tj) = tcM(tj) + 
t-DM{tj){Xi/\-REF)'^ + tm{tj); (3) the variances of the sam- 
ples twi{tj) are known. 
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Figure 1. Average power spectra of pre- and post-correction timing residuals, in ttie 20-cm band, with four combinations of signals. 
The solid line shows the pre-correction spectrum and dashed line shows the post-correction spectrum. For the cases where variations 
in DM are included in the simulation, the pre-correction spectrum without DM variations is shown with a dotted line. Here the fitting 
routine uses the DM(t) interpolated fitting routine, without fitting a common-mode signal. The vertical bars on the left of each panel 
show the 90% spectral estimation uncertainty for a single realisation (left- most bar) and the average of fOOO realisations (right bar). 




Lag (days) 

Figure 4. Cross-correlation of post-correction residuals with 
input model GWB, for the simulations representing PSR 
J1909— 3744. Solid lines show data from the 20-cm wavelength 
and dashed lines show data from the 50cm-cm band. The cor- 
rection was computed with and without fitting for the common 
mode, indicated by thick and thin lines respectively. The auto- 
correlation of the input GWB is plotted as a dotted line, but it 
is completely obscured by the heavy solid line for the cross cor- 
relation in the 20-cm band. 



If Assumption 1 does not hold and we fit for tcm{tj), 
then our method would be sub-optimal. However in any pul- 
sar timing experiment, we must first assume that there is 
a common-mode signal present. If tcmitj) is weak or non- 
existent then we will have a very low corner frequency and 
eflectively we will not fit for tcmitj)- So this assumption is 
tested in every case. 

Assumption 2 will fail if there are wavelength depen- 
dent terms which do not behave like A^, for example the 
scattering effect which behaves more like A*. If these terms 
are present they will corrupt the DM estimate and some 
scattering effects from longer wavelengths may leak into the 
shorter wavelengths due to the correction process. However 
the correction process will not remove any common-mode 
signal, so signals of interest to PTAs will survive the DM 
correction unchanged. 

Assumption 3 will not always be true a priori. Recent 
analysis of single pulses from bright MSPs have shown that 
pulse-to-pulse variations contribute significant white noise 
in excess of that expected from the formal ToA measure- 
ment uncertainty (??). Indeed, many pulsars appear to show 
some form of additional white noise which is currently un- 
explained but could be caused by the pulsar, the interstellar 
medium, or the observing system (see e.g. ??). In any case, 
we cannot safely assume that the uncertainties of the timing 
residuals aij accurately reflect the white noise level. If the 
(Jij are incorrect, our fit parameters, tdt) and tdit), will no 
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Figure 2. As for Figure[T] except the fitting routine uses tlie DM(t) interpolated fitting routine in addition to the wavelength independent 
signal, C(t). 
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Figure 3. Average power spectra of the error in DM(t) after fitting to simulations with realistic sampling and uncertainties. The 
simulations contained white noise, DM variations and, in the left panel, a model GWB. The solid black line shows the power spectrum 
of DMeat(i) — DMin(i). The dotted line is the power spectrum of the white noise only. The dashed line is the post-correction power 
spectrum of the residuals, after subtracting the model GWB signal if present. The crosses mark the sum of the black line and the dotted 
line. 
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Table 2. Scintillation and dispersion properties for the 20 PPTA 
pulsars, at a reference wavelength of 20 cm. The scintillation 
bandwidth (I'o) and time scale (r^iff) are averaged over a large 
number of PPTA observations except for values in parenthesis 
which are are taken from ?. -Diooo is the value of the structure 
function at 1000 days and Tg is the optimal sampling interval for 
tDM(t)- 



Source 


1^0 


Tdifl 


-Diooo 






(MHz) 


(s) 







J0437-4715 


1000 


2486 


1.6 


0.2 


J0613-0200 


1.64 


4500 


0.3 


1 


J0711-6830 


36 


1962 


1.9 


2 


J1022+1001 


65 


2334 


0.14 


2 


J1024-0719 


268 


4180 


6.2 


1 


J1045-4509 


(0.094) 


(119) 


690 


0.25 


J1600-3053 


0.09 


271 


24 


0.5 


J1603-7202 


5 


582 


5.5 


1 


J1643-1224 


0.022 


582 


65 


0.5 


J1713+0747 


24 


2855 


0.31 


1 


J1730-2304 


12.4 


1615 


20 


1 


J1732-5049 


5.4 


1200 


10.0 


1 


J1744-1134 


60 


2070 


1.3 


1 


J1824-2452A 


(0.025) 


(75) 


250 


0.33 


J1857+0943 


5.5 


1464 


0.9 


2 


J1909-3744 


37 


2258 


3.5 


0.33 


J1939+2134 


1.2 


327 


8.9 


0.33 


J2124-3358 


(1170) 


(10705) 


0.4 


2 


J2129-5721 


17.1 


3060 


0.49 


2 


J2145-0750 


195 


3397 


0.15 


2 



longer be minimum variance estimators; however, they will 
remain unbiased. This means that DM estimation will be 
unbiased and the DM correction will not remove any GWB 
(or other common-mode signal) although the correction pro- 
cess may add more white noise than optimal. It should be 
noted that if all the Oij were changed by the same factor our 
DM correction would be unchanged. Fortunately the actual 
white noise is relatively easy to estimate from the obser- 
vations because there are more degrees of freedom in the 
white noise than in the red noise, so in practise we use the 
estimated white noise rather than the formal measurement 
uncertainties oa. 



5 APPLICATION TO PPTA OBSERVATIONS 

We have applied the new DM correction technique to the 
PPTA data set (?). Observations of the PPTA are made 
in three wavelength bands: '10cm' (~ 3100 MHz); '20cm' 
(~ 1400 MHz); and '50cm cm' (~ 700 MHz). The 10-cm and 
20-cm bands have been constant over the entire time span, 
however the long wavelength receiver was switched from a 
centre frequency of 685 MHz to 732 MHz around MJD 55030 
to avoid RFI associated with digital TV transmissions. To 
allow for changes in the intrinsic pulse profile between these 
different wavelength bands, we fit for two arbitrary delays 
between one wavelength band and each of the other bands. 
However we did not allow an arbitrary delay between 685 
and 732 MHz because the pulse shape does not change sig- 
nificantly in that range. 

We began our analysis by using the procedure described 
in Section [3] to compute pilot estimations of DM(t) and 



tcM{t) for each of the 20 pulsars, using a sampling inter- 
val Ts = 0.25 yr. Figure [5] shows the DM(f) derived from the 
above. Our results are consistent with the measurements 
made by ? for the ~ 500 days of overlapping data, which is 
expected since they are derived from the same observations. 



5.1 Determining the sampling interval 

As discussed in Section |3l we can use the diffractive time 
scale Tdiff to predict the magnitude of the DM variations in 
a given pulsar. This value can be computed directly from 
observations, however it is always quite variable on a day to 
day time scale (see Section [6] for discussion), and for a few 
pulsars Tdiff approaches the duration of the observations, so 
it can be hard to measure. Nevertheless we have obtained 
an estimate of the average Tdiff from the dynamic spectra 
for each pulsar, and this is given in Table [2] We have not 
provided an error estimate on the average Tdiff because the 
variation usually exceeds a factor of two so the values tabu- 
lated are very rough estimates. 

We also computed the structure function Dtdm directly 
from the turn values. These structure functions, scaled to 
delay in /xs^ at 20 cm, and those estimated from Tdiff, are 
shown in Figure |S] The value Dtdm(1000 days) is given in 
Table [3 

For each pulsar we also make an estimate of the white 
noise power directly from the power spectrum of the resid- 
uals. The estimates of at each wavelength are given in 
Table [U 

We then use the Dtdm estimates and Equation (|12|l to 
generate a model power spectrum PtdimC/) at the reference 
wavelength (20 cm) for each pulsar. These assume a Kol- 
mogorov spectral exponent. From these model spectra and 
the corresponding Pw,tdm, tabulated in Tabled we deter- 
mine the corner frequency and the corresponding sample 
interval Ts for DM for each pulsar. As we do not have any 
a priori knowledge of tcm for the PPTA pulsars we choose 
the same sample interval for £cm as for iDM- 



5.2 Results 

The measured DM(t) sampled at the optimal interval Ts 
is overlaid on the plot of the pilot analysis with Ts = 
0.25 yr on Figure [S] It is not clear that there are measur- 
able variations in DM in PSRs J1022-f 1001, J2124-3358 or 
J2145— 0750, but one can see that there are statistically sig- 
nificant changes with time for the other pulsars. In general, 
the 'optimally sampled' time series (dashed line) follows the 
DM trend with less scatter. However, there are some sig- 
nificant DM fluctuations that are not well modelled by the 
smoother time-series. In particular we do not model the sig- 
nificant annual variations observed in PSR J0613— 0200, and 
we must add a step change to account for the 250 day in- 
crease in DM observed in PSR J1603— 7202 (these features 
are discussed more fully in Section [Sjl . These variations do 
not follow the Kolmogorov model that was used to derive the 
optimal sampling rate, and therefore we must use a shorter 
Ts so we can track these rapid variations. These results il- 
lustrate the importance of making a pilot analysis before 
deciding on the sample interval. The ISM is an inhomoge- 
neous turbulent process and an individual realisation may 
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Figure 5. DM as a function of time for 20 PPTA pulsars. Solid lines show values measured with intervals of 0.25 yr. In the cases where 
the optimal Ts at a wavelength of 20cm is longer than 0.25 yr, a dashed line is added showing DM(ti) measured with this time step. 
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Figure 6. Structure functions of dispersive delay at 20 cm. The square markers indicate the structure function as measured directly 
from the DM time-series in Figure |5] error bars are derived by simulation of white noise. The solid lines show the extrapolation from 
the scintillation time scale Tdiff assuming a Kolmogorov spectrum, dashed lines mark the region occupied by 66% of simulated data sets 
having Kolmogorov noise with the same amplitude. The dotted lines show a Kolmogorov spectrum with the amplitude set to match the 
real data at a lag of 1000 days. 
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Table 3. Impact of the DM corrections on the timing parame- 
ters, as determined in the 20-cm band. For each pulsar we present 
the change in u and u due to the DM correction, relative to the 
measurement uncertainty, and the ratio of the RMS of the residu- 
als before (Eprc) and after (Spost) DM correction. Also included 
is the ratio of the power spectral density before (Ppre) and after 
(^post) correction, averaged below fc- The final column indicates 
if we believe that the DM corrections have 'improved' the data 
set for the purpose of detecting common- mode red signals. 



PSR 




|Ai> 


Spost 


-Ppost 


Imp. 








Ppro 


J0437-4715 


92 


48 


0.6 


0.15-0.25 


Y 


J0613-0200 


0.16 


2.9 


1.1 


0.3-1.2 


y 


J0711-6830 


3.9 


5.5 


1.0 


0.4-1.6 


y 


J1022-I-1001 


1.4 


0.3 


1.0 


0.6-2.6 


n 


J1024-0719 


1 


0.91 


1.0 


0.2-0.7 


Y 


J1045-4509 


28 


11 


0.7 


0.22-0.39 


Y 


J1600-3053 


35 


0.51 


1.0 


0.4-0.8 


Y 


J1603-7202 


2.4 


2.5 


1.0 


0.2-0.9 


Y 


J1643-1224 


11 


0.73 


1.7 


1.3-3.1 


N 


J1713-I-0747 


3.2 


6.2 


1.0 


0.2-0.7 


Y 


J1730-2304 


6.5 


1.8 


1.1 


0.9-3.2 


n 


J1732-5049 


2.6 


2.8 


1.0 


0.4-1.4 


y 


J1744-1134 


5.4 


0.48 


1.0 


0.5-2.0 


n 


J1824-2452A 


24 


31 


0.7 


0.29-0.56 


Y 


J1857-I-0943 


4.3 


1 


1.0 


0.2-1.0 


y 


J1909-3744 


28 


5 


1.0 


0.44-0.79 


Y 


J1939-I-2134 


13 


1.7 


0.7 


0.34-0.67 


Y 


J2124-3358 


0.25 


0.056 


1.0 


0.5-1.9 


y 


J2129-5721 


3 


2.1 


1.1 


0.7-2.8 


n 


J2145-0750 


0.22 


0.18 


1.0 


0.2-1.0 


y 



not behave much like the statistical average. The DM(t) for 
PSR J1909— 3744 is also instructive. It is remarkably linear 
over the entire observation interval. This linearity would not 
be reflected in the timing residuals at a single wavelength 
because a quadratic polynomial is removed in fitting the 
timing model. It can only be seen by comparing the residu- 
als at different wavelengths. Such linear behaviour implies a 
quadratic structure function and a power spectrum steeper 
than Kolmogorov. 

5.3 Performance of DM Correction 

The simplest and most widely used metric for the quality of 
timing residuals is the RMS of the residuals. Thus a natural 
measure of the performance of DM correction would be the 
ratio of the RMS of the 20-cm residuals before and after DM 
correction. This ratio is provided in Table [S] However, for 
most of these pulsars, the RMS is dominated by the white 
noise and so does not change appreciably after DM correc- 
tion. Furthermore much of the effect of DM(t) variations is 
absorbed by fitting for the pulse frequency and its derivative. 
Thus the ratio of the RMS before and after DM correction 
is not a very sensitive performance measure. As noted by 
?, the DM correction has a significant effect on the pulsar 
spin parameters, which can give an indication of the mag- 
nitude of the DM correction. Table [3] lists the change in v 
and as a factor of the measurement uncertainty, caused 
by applying the DM correction. However, there are system- 
atic uncertainties in the estimation of the intrinsic values of 



V and ii that may be greater than the error induced by DM 
variations. 

Judging the significance of the DM corrections depends 
on the intended use of the data set. Since a major goal of the 
PPTA is to search for common-mode red signals, we choose 
to consider the impact of the DM corrections on the low fre- 
quency noise. In principal, the DM correction should reduce 
the noise at frequencies below /c, and therefore we have 
estimated the ratio of the pre- and post-correction power 
spectrum of the 20-cm residuals, averaged over all frequen- 
cies below /c. We caution that the spectral estimates are 
highly uncertain, and for many pulsars we average very few 
spectral channels so the error is non-Gaussian. Therefore, 
we present these ratios in Table [3] as an estimated 66% un- 
certainty range, determined assuming the spectral estimates 
are x^-distributed with mean and variance equal to the mea- 
sured mean power spectral density. 

There are 9 pulsars for which the DM correction ap- 
pears to significantly reduce the low frequency noise, and 
therefore increases the signal-to-noise ratio for any common- 
mode signal in the data. These pulsars are listed with a 'Y' 
in Table O There are 10 pulsars for which the change in low 
frequency power is smaller than the uncertainty in the spec- 
tral estimation and so it is not clear if the DM correction 
should be performed. Table [3] indicates these pulsars with a 
'y' or 'n', with the former indicating that we believe that the 
DM correction is likely to improve the residuals. However, 
the DM correction fails to 'improve' PSR J1643-1224 under 
any metric, even though we measure considerable DM vari- 
ations (see Figure [5}. As discussed in Section |6] we believe 
that this is due to variations in scattering delay entering the 
DM correction and adding considerable excess noise to the 
corrected residuals. 



6 SCATTERING AND DM CORRECTION 

The most important effect of the ISM on pulsar timing is 
the group delay caused by the dispersive plasma along the 
line of sight. However small scale fluctuations in the ISM 
also cause angular scattering by a diffractive process. This 
scattering causes a time delay to « 0.59qL/c, where 6^0 is 
the RMS of the scattering angle and L is the distance to 
the pulsar. This can be significant, particularly at longer 
wavelengths, because it varies much faster with A than does 
the dispersive delay - approximately as A''. In homogeneous 
turbulence one would expect this parameter to be relatively 
constant with time. If so, the delay can be absorbed into the 
pulsar profile and it will have little effect on pulsar timing. 
However if the turbulence is inhomogeneous the scattering 
delay may vary with time and could become a significant 
noise source for pulsar timing. We can study this effect us- 
ing the PPTA pulsar PSR J1939-I-2134. Afthough this pul- 
sar is unusual in some respects, the scattering is a property 
of the ISM, not the pulsar, and the ISM in the direction 
of PSR J1939-I-2134 can be assumed to be typical of the 
ISM in general. PSR J1939-I-2134 is a very strong source 
and the observing parameters used for the PPTA are well- 
suited to studying its interstellar scattering. The time delay, 
to, can be estimated from the bandwidth of the diffractive 
scintillations, vq, in a dynamic spectrum using the relation- 
ship to — l/2iivo. In fact it is extremely variable, as can be 
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seen in Figure [T] The RMS of to (52 ns at 20 cm) is about 
28% of the mean. We can expect this to increase by a fac- 
tor of (1400 MHz/700 MHz)'' = 16 at 50cm cm. Thus in the 
50cm-cm ToAs there will be delays with RMS variations 
of ~ 830 ns, which do not fit the dispersive A'^ behaviour. 
This will appear in the estimate of toivi at 20 cm, attenuated 
by a factor of ((1400 MHz/700 MHz)^ - 1) = 3 (Equation 
[2|. Therefore the DM correction will bring scattering noise 
from the 50-cm band to the 20-cm band with RMS variation 
~ 270 ns, 5.3 times larger than the scattering noise intrinsic 
to the 20-cm observations. This analysis is corroborated by 
the structure function of DM for this pulsar shown in Fig- 
ure ??, which shows a flattening to about 1 ^s^ at small time 
lags. This implies a white process with RMS variations of 
about 500 ns, consistant with that expected from scattering. 

We have correlated the variations in to with the 20-cm 
residuals before correction and find 18% positive correlation. 
This is consistent with the presence of 52 ns of completely 
correlated noise due to to added to the ToA measurement 
uncertainty of the order of 200 ns. PSR J1939-f 2134 is known 
to show ToA variations that are correlated with the intensity 
scintillations (?) but are much stronger than expected for 
homogeneous turbulence (?). Thus we are confident that the 
observed variation in to is showing up in the 20-cm residuals. 
We expect that contribution to increase in the DM corrected 
residuals to about 300 ns. However this is very difficult to 
measure directly because the DM correction is smoothed 
and the 50cm-cm observations are not simultaneous with 
the 20-cm observations. 

This effect increases rapidly with longer wavelength. If 
we had used 80-cm observations for DM correction, the RMS 
at 80 cm would have been ~ 10 /xs and this would have been 
reduced by a factor of 12 to an RMS of 800 ns in the cor- 
rected 20 cm residuals. Clearly use of low frequency antennas 
such as GMRT (?) or LOFAR (?) for correcting DM fluctu- 
ations in PTAs will have to be limited to weakly scattered 
pulsars. This is an important consideration, but it should 
be noted that the four PPTA pulsars that provide the best 
timing are all scattered much less than J1939-f2134 - all 
could be DM corrected with 80-cm observations or even with 
longer wavelengths. On the other hand there are four PPTA 
pulsars that are scattered 20 to 80 times more strongly than 
J1939-I-2134 and even correction with 50cm-cm data causes 
serious increases in the white noise. 

An extreme example is PSR J1643— 1224. Under the 
above assumption, the expected white noise (2.0 pis) due to 
scattering at 20 cm exceeds the radiometer noise (0.63 /is). 
The white scattering noise at 50cm cm is much larger (32 
/is) and about a third of this makes its way into the DM- 
corrected residuals at 20 cm. This is also corroborated by 
the structure function for this pulsar in Figure ??, which 
shows a flattening to about 10 at small lags. This im- 
plies a white process with RMS variation of ~ 3 /^s which is 
consistant with that expected from scattering. Indeed, this 
pulsar is the only pulsar with significant DM variations for 
which the DM correction increases the noise in the 20-cm 
residuals under all metrics. It is also important to note that 
observing this source at the same frequencies with a more 
sensitive telescope will not improve the signal to noise ra- 
tio, because the noise, both before and after DM-correction, 
is dominated by scattering. However using a more sensitive 
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Figure 7. DifTractive scattering delay, io, measured from scin- 
tillation bandwidth, vq, in observations of PSR J1939+2134 at a 
wavelength of 20 cm. The error bars are derived from the fit for 
and so are roughly proportional to io- 



telescope could improve matters by putting more weight on 
observations at 10 cm, where scattering is negligible. 

Finally however, we note that the usefulness of long 
wavelength observations would be greatly improved if one 
could measure and correct for the variation in scattering de- 
lays. This may be possible using a technique such as cyclic 
spectroscopy, however this has only been done in ideal cir- 
cumstances and with signal-to-noise ratio such that individ- 
ual pulses are detectable (?). It is still unclear if such tech- 
niques can be generalised to other observations, or if this 
can be used to accurately determine the unscattered ToA. 



7 SCHEDULING FOR DM CORRECTION 

If there were no DM variation, one would spend all the ob- 
serving time at the wavelength for which the pulsar has the 
greatest ToA precision (see ? for discussion on choice of ob- 
serving wavelength). The reality is, of course, that we need 
to spend some of the time observing at other wavelengths 
to correct for DM variations. In this section we will present 
a strategy for choosing the observing time at each wave- 
length, attempting to optimise the signal to noise ratio of 
the common-mode signal, tcM ■ We take the PPTA observa- 
tions at Parkes as our example, but this work can easily be 
generalised to any telescope. 

At Parkes it is possible to observe at wavelengths of 10 
and 50cm cm simultaneously because the ratio of the wave- 
lengths is so high that the shorter wavelength feed can be 
located co-axially inside the longer wavelength feed. How- 
ever the 20-cm receiver does not overlap with either 10 or 
50cm cm and so must be operated separately. 

As noted earlier we can write tcM = X^i ^i^oi so the vari- 
ance of tcM is given by o'^Qy^ = '^^hla^. However, the so- 
lution is not linear in the terms because they also appear 
in the bi coefficients. At present, observing time at Parkes is 
roughly equal between the two receivers, so we use the exist- 
ing power spectral densities as the reference. We will assume 
that the total observing time is unity and the time devoted 
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to 10 and 50cm cm, which are observed simultaneously, is 
< « < 1. 

The variance of any toi is inversely proportional to the 
corresponding observing time. We use 20 as the refer- 
ence because it usually has the smallest TOA uncertainty 
in PPTA observations. Therefore, we define ctIq = 1/(1 — a) 
as the reference. Then we assume that, with equal observing 
time tJio = xa2o and aso — y(T2o, so as scheduled we would 
have (Tio = x/a and uso — y/a. We can then determine the 
increase in white noise caused by correcting for dispersion as 
a function of a, the time devoted to 10 and 50cm cm. The re- 
sults are shown for all the PPTA pulsars in Table [T] One can 
see that all the pulsars are different and the optimal strate- 
gies range from a ~ 0.2 to a ~ 1.0 (i.e. 100% of time spent 
using the dual 10 and 50cm cm system). For the four 'best' 
pulsars, PSRs J0437-4715, J1713-h0747, J1744-1134, and 
J1909-3744, the optimal strategy has a > 0.7. 

This suggests that a useful improvement to PTA perfor- 
mance could come from deploying broadband receivers, so 
that correction for DM(t) can be done with a single observa- 
tion. This also has the benefit of reducing the difficulties of 
aligning pulsar profiles measured with different receivers at 
different times, and would therefore allow for more accurate 
measurement of DM variations. 



8 THE INTERSTELLAR MEDIUM 

The PPTA DM{t) observations provide an interesting pic- 
ture of the ionised ISM on au spatial scales. The overall 
picture can be seen in Figures [5] and (6] In Figure [5] it is ap- 
parent that 17 of the 20 PPTA pulsars have measurable 
DM(t) variations. In Figure [6] it can be seen that 13 of 
these 17 show power-law structure functions, as expected 
in the ensemble average Of these, 8 are roughly consistant 
with an extrapolation from the diffractive scales at the Kol- 
mogorov spectral exponent, an average dynamic range of 
4.8 decades. However 5 are considerably higher than is pre- 
dicted by a Kolmogorov extrapolation. They may be locally 
Kolmogorov, i.e. an inner scale may occur somewhere be- 
tween the diffractive scale and the directly measured scales 
of 100 to 2000 days, but establishing this would require a 
detailed analysis of the apparent scintillation velocity which 
is beyond the scope of this paper. Two of these five pulsars, 
J1045— 4509 and J1909— 3744, were already known to be in- 
consistent with a Kolmogorov spectral exponent (?), and it 
is clear, with the additional data that are now available, that 
J1024-0719, J1643-1224 and J1730-2304 should be added 
to this list. When the spatial power spectrum of a stochas- 
tic process is steeper than the Kolmogorov power-law, it can 
be expected to be dominated by linear gradients and show 
an almost quadratic structure function. Indeed inspection of 
Figure [5] shows that the 5 steep spectrum pulsars all show a 
strong linear gradient in DM{t). 

It is interesting that two of the five steep spectrum pul- 
sars, J1730— 2304 and J1909— 3744, show structure functions 
which are clearly steeper than Kolmogorov in the observed 
range and are thus converging towards a Kolmogorov spec- 
trum at smaller time lags. The three other steep spectrum 
pulsars do not show this behavior in the structure functions 
for two reasons: (1) J1024— 0719 appears to flatten, appar- 
ently due to a white noise contribution which is comparable 



with the error bars and is probably due to underestimation 
of bias correction due to the errors; (2) J1045— 4509 and 
J1643— 1224 are highly scattered and are showing the effect 
of a scattering contribution at small lags. The scattering 
contribution to J1730-2304 and J1909-3744 is negligible. 

The time series DM(t) shown in Figure [5] often show 
behaviour that does not look like a homogeneous stochastic 
process. For example, PSR J1603— 7202 shows a large in- 
crease for ~250 days around MJD 54000 and J0613-0200 
shows clear annual modulation. The increase in DM for 
J1603— 7202 suggests that a blob of plasma moved through 
the line of sight. If we assume the blob is halfway between 
the pulsar and the Earth, the line of sight would have moved 
by about 0.5 au in this time, and if the blob were spherical it 
would need a density of ~200 cm^^. This value is high, but 
comparable to other density estimates for au-scale structure 
based on 'extreme scattering events' (??). 

We computed the power spectra of DM{t) for all the 
pulsars to see if the annual modulation that is clear by 
eye in PSR J0613— 0200 is present in any of the other pul- 
sars. For four pulsars we find a significant (> 5-cr) detection 
of an annual periodicity, PSRs J0613-0200, J1045-4509, 
J1643-1224 and J1939+2134. 

The most likely explanation for the annual variation in 
T)M{t) is the annual shift in the line of sight to the pulsar 
resulting from the orbital motion of the Earth. The trajec- 
tory of the line of sight to three example PPTA pulsars are 
shown in Figure [S] The relatively low proper motion and 
large parallax of the PPTA pulsars means that the trajec- 
tory of the line of sight to many of the PPTA pulsars show 
pronounced ripples. However, unless the trajectory is a tight 
spiral, the annual modulation will only be significant if there 
is a persistent gradient in the diffractive phase screen. 

The presence of persistent phase gradients and annual 
modulation in J1045-4509 and J1643-1224 is not surpris- 
ing because the ISM associated with each of these pulsars 
has a steeper than Kolmogorov power spectrum. Indeed, 
the measured DM(t) for these pulsars do show a very lin- 
ear trend, which in itself evidence for a persistent phase 
gradient. The other steep spectrum pulsars, J1024— 0719, 
J1730-2304 and J1909-3744, have higher proper motion, 
which reduces the amplitude of the annual modulation rela- 
tive to the long term trend in DM(t). We note that the spec- 
tral analyses for PSRs J1024-0719 and J1909-3744 suggest 
annual periodicities, and it may be possible to make a sig- 
nificant detection by combining the PPTA data with other 
data sets. 

PSR J1939-I-2134 does not show a steep spectrum, how- 
ever its proper motion is very low compared to its parallax, 
and therefore the trajectory spirals through the ISM, reduc- 
ing the requirement for a smooth phase screen. The annual 
modulation of J0613— 0200 may be somewhat different, since 
it does not have a steep spectrum and although the proper 
motion is small the trajectory does not spiral (see Figure [S}. 
This suggests that for J0613— 0200 the turbulence could be 
anisotropic with the slope of the gradient aligned with the 
direction of the proper motion. Anisotropic structures are 
believed to be quite common in the ISM (??). However one 
can imagine various other ways in which this could occur, 
particularly in an inhomogeneous random process, and inho- 
mogeneous turbulence on an au spatial scale is also believed 
to be common in the ISM (???). 
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Persistent spatial gradients will cause a refractive shift 
in the apparent position of the pulsar, and because of dis- 
persion the refraction angle will be wavelength dependent. 
This refractive shift appears in the timing residuals as an 
annual sine wave which changes in amplitude like A'^. When 
the DM(t) is corrected this sine wave disappears and the in- 
ferred position becomes the same at all wavelengths. These 
position shifts are of order 10"" (A/20 cm) arcseconds for 
all four pulsars. 

Note that the trajectory of the lines of sight shown on 
Figure [8] may appear quite non-sinusoidal, but the annual 
modulation caused by the Earth's orbital motion in a linear 
phase gradient will be exactly a sine wave superimposed on a 
linear slope due to proper motion. This will not generate any 
higher harmonics unless the structure shows significant non- 
linearity on an au scale. We do not see second harmonics of 
the annual period, which suggests that the spatial structure 
must be quite linear on an au scale. 

Annual variations in DM are also observed in pulsars 
for which the line of sight passes close to the Sun because 
of free electrons in the solar wind (??). In the PPTA, a 
simple symmetric model of the solar wind is used to remove 
this effect, but this is negligible for most pulsars. For the 
three pulsars where it is not negligible, the effect of the solar 
wind persists only for a few days at the time when the line 
of sight passes closest to the Sun. Neither the magnitude, 
phase nor shape of the variations seen in our sample can 
be explained by an error in the model of the solar wind. 
Changes in ionospheric free electron content can be ruled 
out for similar reasons. 

In summary the ISM observations are, roughly speak- 
ing, consistent with our present understanding of the ISM. 
However the data will clearly support a more detailed anal- 
ysis, including spectral modelling over a time scale range in 
excess of 10^ from the diffractive scale to the duration of the 
observations. It may also be possible to make a 2 dimensional 
spatial model of the electron density variations for some of 
the 20 PPTA pulsars, although such detailed modelling is 
far beyond the scope of this paper. Preliminary attempts 
to model the DM variations in PSR J0613— 0200, assuming 
that the DM can be approximated as a linear gradient dur- 
ing the observation period, suggest that such modelling may 
be useful for both studying the ISM and to improve the DM 
correction. 



9 CONCLUSIONS 

We find that it is necessary to approach the problem of es- 
timating and correcting for DM(f) variations iteratively, be- 
ginning with a pilot analysis for each pulsar and refining that 
analysis as the properties of that pulsar and the associated 
ISM become clearer. Each pulsar is different and the ISM in 
the line of sight to each pulsar is different. The optimal anal- 
ysis must be tailored to the conditions appropriate for each 
pulsar and according to the application under consideration. 

We sample the DM(t) just often enough that the varia- 
tions in DM are captured with the minimum amount of addi- 
tional white noise. Likewise, we must also sample a common- 
mode signal tcM{t) at the appropriate rate. In this way we 
can correct for the DM variations at frequencies where it is 
necessary, and we can include tcM(t) at frequencies where 
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Figure 8. Trajectories through the ISM of the line of sight to 
PSRs J0613-0200 (dashed Hne), J1643-1224 (soUd black line) 
and J1909— 3744 (grey line). It was assumed that the scattering 
takes place half way between the pulsar and the Earth and the 
motion of the plasma was neglected. The trajectories are marked 
with a cross at the DM sampling interval of 0.25 yr. 

it is necessary, but not fit for either at frequencies where the 
signal is dominated by white noise. 

By including the common-mode signal in the analysis 
we preserve the wavelength-independent signals of interest 
for pulsar timing arrays and we improve the estimate of 
the pulsar period and period derivative. Without estimat- 
ing the common mode, a significant fraction of wavelength- 
independent signals, such as: errors in the terrestrial clocks; 
errors in the planetary ephemeris; and the effects of gravita- 
tional waves from cosmic sources, would have been absorbed 
into the DM correction and lost. 

We have applied this technique to the PPTA data set, 
which improves its sensitivity for the detection of low fre- 
quency signals significantly. The estimated DM{t) also pro- 
vides an unparallelled measure of the au scale structure of 
the interstellar plasma. In particular it confirms earlier sug- 
gestions that the fluctuations often have a steeper than Kol- 
mogorov spectrum, which implies that an improved physical 
understanding of the turbulence will be necessary. We also 
find that persistent phase gradients over au scales are rel- 
atively common and are large enough to cause significant 
errors in the apparent positions of pulsars unless DM cor- 
rections are applied. 
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APPENDIX A: CONSTRAINED LEAST 
SQUARES FITTING IN TEMP02 

The least squares problem of fitting tlie timing model to the 
residuals can be written in matrix form as 

R = MP + E. (Af) 

Here R is a column vector of the timing residuals, P is 
a column vector of fit parameters, including DM(tj) and 
tcM{tj) as well as the other timing model parameters. M 
is a matrix describing the timing model and E is a column 
vector of errors. The least-squares algorithm solves for P, 
matching MP to R with a typical accuracy of E. 

The sampled time series DM(fj) and tcMitj) are co- 
variant with the timing model, so they must be con- 
strained to eliminate that covariance or the least squares 
solution will fail to converge on a unique solution. These 
constraints have the form of linear equations of DM(tj) 
and tcmitj), such as: Y^DM{tj) = 0; I]tcM(tj) = 0; 
J2tjtcM{tj) = 0, E^iteMfe) = 0; EsinMj)tcMfe) = 0; 
'^cos{ujtj)tcM{tj) ~ 0; etc. Augmented with these equa- 
tions, the least-squares problem becomes 
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where B is a matrix describing the constraints, e is a column 
vector of weights for the constraints. In our case C = 0, 
though it need not be in general. The least-squares solution 
will then find a vector P that matches both MP to R, with 
a typical accuracy of E, and also matches BP to C, with a 
typical accuracy of e. By making e very small we can enforce 
the constraints with high accuracy. This scheme has been 
called 'the method of weights' (?). 

If the uncertainties in the estimates of DM(tj) and 
tcmitj) are not expected to be equal, for instance if the 
different observing wavelengths are irregularly sampled and 
the ToA uncertainties are variable across sampling windows, 
then it can be advantageous to use weighted constraints. 
Then the constraints take the form ^ WjDM{tj) = 0, and 
we need to estimate the uncertainties of the parameters to 
obtain the optimal weights. These uncertainties can be de- 
termined from the least-squares solution in which the timing 
residuals are described purely by Equation Q. This problem 
is linear and the covariance matrix of the parameters can be 
written in closed form without even solving for the parame- 
ters. The diagonal elements of the covariance matrix are the 
variances of the parameters and the weights, Wj, are the 
inverse of the square roots of the corresponding variances. 



If y(t) is real valued then by the Wiener-Khinchin theorem, 

C(r)= / cos(2^/r)P(/)d/, (B3) 
Jo 

where P{f) is the one-sided power spectral density of y{t). 
Thus we can then write the structure function in terms of 
the power spectral density as 

D{t) =A 2(1 - cos(27r/r))P(/) d/, (B4) 
Jo 

It should be noted that this expression for D{t) is valid if 
D{t) exists. It is not necessary that C(r) exist. For the case 
of a power-law, P{f) = Af~°', we can change variables using 
X = /r, and obtain 

roo 

D(t)=t°'-^A 2(l-cos{2-Kx))x~°'dx. (B5) 
Jo 

The integral (Int) above converges if 1 < q < 3, yielding 

Int. = 2"tv°'-^ sin(-a7r/2)r(l - a), (B6) 

where F is the Gamma function. Thus for Kolmogorov tur- 
bulence, with exponent a — 8/3, we have Int ~ 89.344 and 
the power spectrum can be written 

P(/) ~ 0.0112D(r)r~^/V"*''^- (B7) 



APPENDIX B: RELATION BETWEEN THE 
STRUCTURE FUNCTION AND 
POWER-SPECTRAL DENSITY 

The structure function D{t), of a time series y{t), is well 
defined if y{t) has stationary differences 

D(T) = {{y{t)-y{t + T)f). (Bl) 

If y{t) is wide-sense stationary D{t) can be written in terms 
of the auto covariance C(r) by expansion of the square 

D{t) = {y{tf) + {y{t + Tf)-2{y{t)y{t + T)) 

= 2(C(0)-C(r)) (B2) 



